2012/10/17

PDL, PDL::Stats 예제 t분포

t-분포정규분포와 유사한 확률분포이며, 모집단의 표준편차를 알지 못하며(표본의 표준편차로 모집단의 표준편차를 대용) 표본의 크기가 작을 때(n<30) 주로 이용한다. 자유도(n-1)가 낮을수록 봉이 낮아지고 꼬리가 높아지며 자유도가 높을수록 정규분포에 가까워진다.

#!/usr/bin/env perl
use 5.010;
use PDL;
use PDL::Stats;
use PDL::GSL::CDF;
use PDL::Graphics::PGPLOT::Window;
use Math::GSL::Randist qw/gsl_ran_tdist_pdf/;

# t 분포
my $win = pgwin();
my $x = zeroes(100)->xlinvals(-4, 4);
$win->env(-4,4,0,1,{Title=>'t distribution'});
$win->hold();
#CDF
$win->line($x, gsl_cdf_tdist_P($x, 29),{COLOR=>'RED'});
$win->line($x, gsl_cdf_tdist_P($x, 4),{COLOR=>'MAGENTA'});
#PDF
$win->line($x, [ map { gsl_ran_tdist_pdf($_, 29) } $x->list ],{COLOR=>'BLUE'});
$win->line($x, [ map { gsl_ran_tdist_pdf($_, 4) } $x->list ],{COLOR=>'CYAN'});

# t분포 자유도10, 우측꼬리확률 0.05일때 t확률변수 값은?
say gsl_cdf_tdist_Pinv(1-0.05, 10);
# t분포 자유도10, 양측꼬리확률 0.05일때 t확률변수 값은?
say gsl_cdf_tdist_Pinv(pdl(0.05/2, 1-0.05/2), 10);


1.81246112281168
[-2.2281389  2.2281389]

자유도29(표본수 30), 자유도4(표본수 5)에 대해서 PDF그래프를 보면 자유도가 낮을수록 펑퍼짐하며 높을수록 봉이 높아지며  정규분포에 가까워짐을 볼 수 있다. 분포가 펑퍼짐하다는 것은 표본수가 작을 수록 불확실성이 크기때문에 양쪽으로 치우칠 확률이 더 높다는 것을 의미한다.

2012/10/12

PDL, PDL::Stats 예제 정규분포

정규분포(가우스분포)는 앞서 올린 이항분포의 n이 아주 커질 경우 근사적으로 따르게되는 확률분포이며(참고) 통계학에서 기본이 되는 연속확률변수 이다. 평균과 표준편차가 다른 많은 정규분포가 존재할 수 있으나 z점수를 통해 평균 0, 표준편차 1 인 표준정규분포로 변환해서 계산가능하다. 계산기 및 컴퓨터가 발달하지 않았던 옛날에는 주로 미리 계산된 표를 통해서 계산했지만 요즘에서는 굳이 그럴 필요없이 컴퓨터에 미리정의된 함수를 활용하면 간단하다.

다음은 표준정규분포의 확률밀도함수(PDF)누적분포함수(CDF)를 그리고 정규분포을 이용해서 몇가지 예에 대해 계산하는 코드이다.

#!/usr/bin/env perl
use 5.010;
use PDL;
use PDL::Stats;
use PDL::GSL::CDF;
use PDL::Graphics::PGPLOT::Window;
use Math::GSL::Randist qw/gsl_ran_gaussian_pdf/;

# 정규분포
my $win = pgwin();
my $x = zeroes(100)->xlinvals(-4, 4);
$win->env(-4,4,0,1,{Title=>'standard normal(gaussian) distribution'});
$win->hold();
#CDF
$win->line($x, gsl_cdf_gaussian_P($x,1), {COLOR=>'RED'});
#PDF
$win->line($x, [ map { gsl_ran_gaussian_pdf($_, 1) } $x->list ], {COLOR=>'BLUE'});
# PDL::Stats::Distr의 pdf_gaussian 함수 사용시
#$win->line($x, $x->pdf_gaussian(0, 1), {COLOR=>'BLUE'});

# 시험결과 학생들의 성적이 평균이 63점 분산이 100인 정규분포를 따른다고 하면
# 50점 미만의 학생은 몇 %인가 ?
# Z = (50-63)/sqrt(100) 이므로
say gsl_cdf_gaussian_P((50-63)/sqrt(100), 1)*100, " %";

# 상위 10%에게 A를 주면 A를 받기 위해서 몇점 이상이 되어야 하는가?
my $Z = gsl_cdf_gaussian_Pinv(0.9, 1); # 표준정규분포에서 하위(100-10)=90% 위치의 Z값
# Z = (X-63)/sqrt(100) 이므로 점수X는
say $Z*sqrt(100)+63


<결과>
9.68004845856103 %
75.815515655446


그래프에서 파란색이 PDF, 빨간색이 CDF의 그래프이다. 표준정규분포의 확률변수구간은 음의무한대에서 양의 무한대 까지이나 -4~4구간만 그렸으며 PDF아래의 면적이1 CDF는 PDF를 구간별로 적분한 값으로 0에서 1에 수렴한다.

2012/10/11

PDL, PDL::Stats 예제 이항분포

앞서 올린 PDL, PDL::Stats 로 기본 통계량 계산하기 에 이어 이항분포에 대한 예를 들어보면 다음과 같다.

이항분포는  연속된 n번의 독립적 시행에서 각 시행이 확률 p를 가질 때의 확률분포를 말하는데 말로 하는 것 보다 다음 그림을 보면 그 의미가 명확해진다.

위 틀의 꼭대기에서 구슬을 떨어뜨리면 각 층에서 좌우로 갈 확률p는 반반(1/2)으로 각층을 통과 하며(독립시행) 구슬이 제각각의 칸으로 떨어지게 되는데 아래 구슬이 쌓인 모양이 바로 이항분포가 되는 것이다.

 그러면 n=30, p=0.5 일때 이항분포의 확률질량함수 PMF[Probability Mass Function](연속확률변수의 PDF에 해당), 누적분포함수 CDF[Cumulative Distribution Function] 의 예를 코딩해 보면 다음과 같다.

#!/usr/bin/env perl
use 5.010;
use PDL;
use PDL::Stats;
use PDL::GSL::CDF;
use Math::GSL::Randist qw/gsl_ran_binomial_pdf/;
use PDL::Graphics::PGPLOT::Window;

my $win1 = pgwin();
my $x = sequence(30);
$win1->env(0,30,0,1,{Title=>'binomial distribution'});
$win1->hold();

#CDF
# gsl_cdf_binomial_P는 PDL::Stats에 포함된 PDL::GSL::CDF 모듈의 함수
$win1->bin($x, gsl_cdf_binomial_P($x, 0.5, 30),{COLOR=>'RED'});

#PMF
# Math::GSL::Randist의 gsl_ran_binomial_pdf 함수를 사용
# PDL용 함수가 아니라서 piddle을 계산하지 못하므로 map으로 piddle $x내부의 값을
# map 함수를 사용하여 리스트를 순회하며 적용
$win1->bin($x, [ map { gsl_ran_binomial_pdf($_, 0.5, 30) } $x->list ], {COLOR=>'BLUE'});
# PDL::Stats::Distr의 pmf_binomial 함수를 사용할 시
# $win1->bin($x, $x->pmf_binomial(30, 0.5), {COLOR=>'BLUE'}); 

#동전을 30번 던졌을때 12번 앞면이 나올 확률
# PMF 그래프 가로축 12일때 세로축의 값
say gsl_ran_binomial_pdf(12, 0.5, 30) * 100, " %";
#say piddle(12)->pmf_binomial(30, 0.5) * 100, " %"; # pmf_binmoial 함수를 사용할 시
#동전을 30번 던졌을때 12번 이하로 앞면이 나올 확률
# CDF 그래프 가로축 12일때 세로축의 값
# 결론적으로 12까지의 PMF그래프 값을 누적한(적분) 것
say gsl_cdf_binomial_P(12, 0.5, 30) * 100, " %";


<결과>


8.05530929937956 %
18.0797304026779 %

위 그래프에서 파란색이 PMF이고 빨간색이 CDF이다. PMF그래프 아래의 면적을 구하면 1이 되고 CDF는 0에서 시작하여 1까지 증가한다는 사실을 보면 확률의 분포와 누적확률이 시각적으로 느껴질 것이라고 본다.

2012/10/10

PDL, PDL::Stats 로 기본 통계량 계산하기

Perl에는 Matlab같은 PDL이란 것이 있다.

여기에 PDL::Stats 모듈을 같이 사용하면 오픈소스 통계패키지인 R같이 여러가지 통계계산도 할 수 있다.

다음은 어떤 데이터의 평균, 분산, 표준편차, 어떤 두 데이터간의 공분산상관계수를 구하는 예제 코드이다.

#!/usr/bin/env perl
use 5.010;
use PDL;
use PDL::Stats;
use PDL::Graphics::PGPLOT::Window;

my $data = pdl( [1,1,5,6,7,3,3,9,4,10,8,6,6,2,6,15] );
say "data: ",$data;
say "average: ",$data->average;
say "variance: ",$data->var;
say "standard deviation: ",$data->stdv;

my $x = pdl( [1,2,3,4,5,6] );
my $y = pdl( [1,1.5,2,2.25,2.5,3] );
my $y2 = pdl( [6,5,4,3,2,1] );

say "x:", $x;
say "y:", $y;
say "y2:", $y2;
say "covariance x,y: ",$x->cov($y);
say "correlation x,y: ",$x->corr($y);
                  
say "covariance x,y2: ",$x->cov($y2);
say "correlation: x,y2: ",$x->corr($y2);

my $win = pgwin();
$win->hold();
$win->env(0,10,0,10);
$win->points($x,$y, {SYMBOL=>'SQUARE', COLOR=>'BLUE'});
$win->points($x,$y2, {SYMBOL=>'PLUS', COLOR=>'CYAN'});

<결과>

data: [1 1 5 6 7 3 3 9 4 10 8 6 6 2 6 15]
average: 5.75
variance: 12.4375
standard deviation: 3.52668399491647
x:[1 2 3 4 5 6]
y:[1 1.5 2 2.25 2.5 3]
y2:[6 5 4 3 2 1]
covariance x,y: 1.10416666666667
correlation x,y: 0.991332701357794
covariance x,y2: -2.91666666666667
correlation: x,y2: -1

















차후에는 여러가지 확률분포에 대해 테스트해보고자 한다.