2.6. 정규분포, 로그정규분포, 스튜던트t분포

2.6.1. 정규분포

정규분포(normal distribution)은 혹은 가우스 정규분포(Gaussian normal distribution)는 값의 제한이 없이 음의 무한대부터 양의 무한대까지 어떤 수든 발생할 수 있지만 특정한 값 주변의 값이 가장 많이 나오는 확률변수의 분포다. 정규분포를 따르는 확률변수의 값은 다음과 같이 대칭인 종(bell) 모양을 이룬다.

2.6.1.1. 정규분포의 확률분포함수

확률변수 \(x\)가 정규분포를 따르는 경우 다음과 같이 표시한다.

\[ x \sim \mathcal{N}(x) \]

위 식에서 \(\mathcal{N}(x)\)은 정규분포의 확률밀도함수를 가리키는 기호다.

정규분포의 확률밀도함수 \(\mathcal{N}(x)\)의 수식은 다음과 같다.

\[ \mathcal{N}(x; \mu, \sigma^2) = \dfrac{1}{\sqrt{2\pi\sigma^2}} \exp \left(-\dfrac{(x-\mu)^2}{2\sigma^2}\right) \]

정규분포의 확률밀도함수는 평균 \(\mu\)와 분산 \(\sigma^2\)이라는 2개의 모수를 가진다. 평균 \(\mu\)는 가장 발생확률이 높은 값으로 정규분포 확률밀도함수의 종 모양의 위치를 지정하는 값이다. 분산 \(\sigma^2\)는 확률밀도함수의 종 모양의 폭을 지정하는 값으로 값이 클수록 폭이 넓어진다. 분산의 역수 \(\sigma^{-2}\)를 정밀도(precision)라고 한다.

다음은 여러가지 모수값을 가지는 정규분포의 확률밀도함수 모양을 나타낸 것이다.

2.6.1.2. 정규분포의 시뮬레이션

scipy 패키지의 stats 서브패키지는 정규분포를 시뮬레이션할 수 있는 norm 클래스를 제공한다. 우선 다음과 같이 모수를 지정한 클래스 인스턴스를 만든다. 이때 loc 인수로 모수 \(\mu\)를 설정하고 scale 인수로 모수 \(\sigma^2\)의 제곱근값 \(\sigma\)를 설정한다.

from scipy.stats import norm

mu = 0
sigma2 = 1
x = norm(loc=mu, scale=np.sqrt(sigma2))
x

<scipy.stats._distn_infrastructure.rv_frozen at 0x13874b2e0>

pdf 메서드를 사용하면 확률밀도함수를 계산할 수 있다. 인수로는 확률변수가 가질 수 있는 값을 넣는다.

xx = np.linspace(-5, 5, 100)
plt.plot(xx, x.pdf(xx))
plt.title("정규분포의 확률밀도함수")
plt.show()

다른 확률분포와 마찬가지로 rvs 메서드를 사용하여 표본측정 시뮬레이션을 할 수 있다.

import numpy as np

np.random.seed(2)

samples = x.rvs(100)
samples

array([-4.16757847e-01, -5.62668272e-02, -2.13619610e+00,  1.64027081e+00,
       -1.79343559e+00, -8.41747366e-01,  5.02881417e-01, -1.24528809e+00,
       -1.05795222e+00, -9.09007615e-01,  5.51454045e-01,  2.29220801e+00,
        4.15393930e-02, -1.11792545e+00,  5.39058321e-01, -5.96159700e-01,
       -1.91304965e-02,  1.17500122e+00, -7.47870949e-01,  9.02525097e-03,
       -8.78107893e-01, -1.56434170e-01,  2.56570452e-01, -9.88779049e-01,
       -3.38821966e-01, -2.36184031e-01, -6.37655012e-01, -1.18761229e+00,
       -1.42121723e+00, -1.53495196e-01, -2.69056960e-01,  2.23136679e+00,
       -2.43476758e+00,  1.12726505e-01,  3.70444537e-01,  1.35963386e+00,
        5.01857207e-01, -8.44213704e-01,  9.76147160e-06,  5.42352572e-01,
       -3.13508197e-01,  7.71011738e-01, -1.86809065e+00,  1.73118467e+00,
        1.46767801e+00, -3.35677339e-01,  6.11340780e-01,  4.79705919e-02,
       -8.29135289e-01,  8.77102184e-02,  1.00036589e+00, -3.81092518e-01,
       -3.75669423e-01, -7.44707629e-02,  4.33496330e-01,  1.27837923e+00,
       -6.34679305e-01,  5.08396243e-01,  2.16116006e-01, -1.85861239e+00,
       -4.19316482e-01, -1.32328898e-01, -3.95702397e-02,  3.26003433e-01,
       -2.04032305e+00,  4.62555231e-02, -6.77675577e-01, -1.43943903e+00,
        5.24296430e-01,  7.35279576e-01, -6.53250268e-01,  8.42456282e-01,
       -3.81516482e-01,  6.64890091e-02, -1.09873895e+00,  1.58448706e+00,
       -2.65944946e+00, -9.14526229e-02,  6.95119605e-01, -2.03346655e+00,
       -1.89469265e-01, -7.72186654e-02,  8.24703005e-01,  1.24821292e+00,
       -4.03892269e-01, -1.38451867e+00,  1.36723542e+00,  1.21788563e+00,
       -4.62005348e-01,  3.50888494e-01,  3.81866234e-01,  5.66275441e-01,
        2.04207979e-01,  1.40669624e+00, -1.73795950e+00,  1.04082395e+00,
        3.80471970e-01, -2.17135269e-01,  1.17353150e+00, -2.34360319e+00])

시뮬레이션된 표본값의 히스토그램을 그리면 다음과 같다.

import seaborn as sns

sns.displot(samples)
plt.xlim(-5, 5)
plt.show()

2.6.1.3. 정규분포의 예

사실 현실세계에서는 상당히 많은 데이터가 정규분포를 따른다. 다음 데이터는 setosa, versicolor, virginica라는 3가지 붓꽃종에 속하는 꽃송이에 대해 꽃잎(petal), 꽃받침(sepal)의 길이(length)와 폭(width)을 0.1cm 단위로 측정한 데이터다.

import statsmodels.api as sm

iris = sm.datasets.get_rdataset("iris").data
iris

	Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species
0	5.1	3.5	1.4	0.2	setosa
1	4.9	3.0	1.4	0.2	setosa
2	4.7	3.2	1.3	0.2	setosa
3	4.6	3.1	1.5	0.2	setosa
4	5.0	3.6	1.4	0.2	setosa
...	...	...	...	...	...
145	6.7	3.0	5.2	2.3	virginica
146	6.3	2.5	5.0	1.9	virginica
147	6.5	3.0	5.2	2.0	virginica
148	6.2	3.4	5.4	2.3	virginica
149	5.9	3.0	5.1	1.8	virginica

150 rows × 5 columns

이 중에서 setosa 종의 꽃받침 길이만 setosa_sepal_length라는 이름으로 뽑아서 분포를 그려보면 정규분포와 유사한 모양을 가지는 것을 볼 수 있다.

setosa_sepal_length = iris[iris.Species == "setosa"]["Sepal.Length"].values

sns.distplot(setosa_sepal_length)
plt.title("붓꽃 꽃받침의 길이 분포")
plt.show()

다음은 팁 데이터에서 지불 금액과 팁 금액의 비율(tip/total_bill)의 분포를 구하면 마찬가지로 정규분포와 유사한 모양을 가진다. 다만 상식적으로는 생각하기 어려운 40% 혹은 70% 근처의 값들이 있는 것을 알 수 있다. 이렇게 다른 값들과 동떨어진 값들을 아웃라이어(outlier)라고 한다.

tips = sns.load_dataset("tips")

tip_ratio = tips["tip"] / tips["total_bill"]

sns.distplot(tip_ratio, rug=True)
plt.title("팁 비율의 분포")
plt.show()

2.6.1.4. QQ플롯과 정규성 검정

그런데 사실 위와 같이 히스토그램 등의 시각적 분포만으로는 데이터의 분포가 정규분포인지 아니면 이와 유사한 다른 분포인지 알기 어렵다. 이를 판별하기 위한 방법으로 두가지를 소개한다.

1. QQ플롯

2. 정규성 검정

QQ플롯은 표본 데이터의 분포와 진짜 정규분포의 분포 형태를 비교하여 표본 데이터가 정규분포를 따르는지 검사하는 간단한 시각적 도구다. QQ플롯을 그렸을 때 직선 모양에 가까울수록 정규분포일 가능성이 높다.

scipy 패키지의 stats 서브패키지는 Q-Q 플롯을 계산하고 그리기 위한 probplot 함수를 제공한다. probplot 함수는 기본적으로 인수로 보낸 데이터 표본에 대한 QQ플롯 정보만을 반환하고 실제 챠트는 그리지 않는다. 만약 차트를 그리고 싶다면 plot 인수에 다음과 같이 matplotlib.pylab 모듈 객체를 넘겨주어야 한다.

다음은 붓꽃 꽃받침 길이에 대해 QQ플롯을 그린 것이다. 직선에 가까운 모양을 띄고 있음을 알 수 있다. QQ플롯의 모양이 계단모양이 된 것은 길이를 측정할 때 0.1cm 단위로 절사했기 때문이다.

sp.stats.probplot(setosa_sepal_length, plot=plt)
plt.show()

팁 비율 데이터에 대해 QQ플롯을 그리면 전반적으로 직선 모양을 따르지만 아웃라이어가 있어서 우측 끝단에서 상단으로 휘어지면서 직선 모양이 깨지는 것을 볼 수 있다.

sp.stats.probplot(tip_ratio, plot=plt)
plt.show()

만약 40%가 넘는 아웃라이어들을 제거하면 직선에 가까운 모양을 띄게 된다.

tip_ratio_ex_outlier = tip_ratio[tip_ratio < 0.4]

sp.stats.probplot(tip_ratio_ex_outlier, plot=plt)
plt.show()

하지만 QQ플롯도 정규분포인지 아닌지 정량적으로 판단해주지는 못한다. 이때는 앞에서 설명한 검정(test)을 사용해야 한다. 어떤 데이터의 분포가 정규분포인지 아닌지를 판단하는 검정을 정규성 검정(normality test)이라고 한다.

정규성 검정을 위한 함수는 scipy 패키지와 statsmodels 패키지에서 제공한다. 정규성 검정 방법이 여러가지이므로 함수도 다양한 함수가 제공한다.

• scipy 패키지가 제공하는 정규성 검정 함수

  • 콜모고로프-스미르노프 검정(Kolmogorov-Smirnov test) : scipy.stats.ks_2samp

  • 샤피로-윌크 검정(Shapiro–Wilk test) : scipy.stats.shapiro

  • 앤더스-달링 검정(Anderson–Darling test) : scipy.stats.anderson

  • 다고스티노 K-제곱 검정(D’Agostino’s K-squared test) : scipy.stats.mstats.normaltest

• statsmodels 패키지가 제공하는 정규성 검정 명령어

  • 콜모고로프-스미르노프 검정(Kolmogorov-Smirnov test) : statsmodels.stats.diagnostic.kstest_normal

  • 옴니버스 검정(Omnibus Normality test) : statsmodels.stats.stattools.omni_normtest

  • 자크-베라 검정(Jarque–Bera test) : statsmodels.stats.stattools.jarque_bera

  • 릴리포스 검정(Lilliefors test) : statsmodels.stats.diagnostic.lillifors

예를 들어 statsmodels 패키지가 제공하는 옴니버스 검정 함수를 사용하여 붓꽃 꽃받침 길이의 정규성 검정을 실시하면 다음과 같이 90.7%의 유의확률을 얻을 수 있다. 정규성 검정의 귀무가설은 “표본 데이터가 정규분포를 따른다”이고 유의확률의 값이 보통 사용하는 1% 혹은 5%의 유의수준보다 높기 때문에 이 귀무가설을 기각할 수 없어 채택하게 된다. 즉, 붓꽃 꽃받침 길이는 정규분포를 따른다.

from statsmodels.stats.stattools import omni_normtest

omni_normtest(setosa_sepal_length)

NormaltestResult(statistic=0.19416346073876928, pvalue=0.907481834276742)

같은 방법으로 아웃라이어를 제거한 팁 비율도 정규분포를 따른다고 판단할 수 있다.

omni_normtest(tip_ratio_ex_outlier)

NormaltestResult(statistic=4.684342675054582, pvalue=0.09611870532870112)

하지만 아웃라이어를 포함하면 유의확률이 0에 가까운 값이 되므로 아웃라이어를 포함한 전체 팁 비율은 정규분포를 따른다고 할 수 없다.

omni_normtest(tip_ratio)

NormaltestResult(statistic=220.8879728815828, pvalue=1.0833932598440914e-48)

2.6.1.5. 정규분포의 모수추정

정규성 검정을 통해 우리가 관심있는 데이터의 분포가 정규분포라는 것을 확인한 후에는 정규분포의 모수 \(\mu\)와 \(\sigma^2\)의 값을 추정해야 한다. 나중에 설명할 최대가능도추정법에 따르면 모수 \(\mu\)의 추정치 \(\hat\mu\)는 표본 데이터의 평균이고 모수 \(\sigma^2\)의 추정치 \(\hat{\sigma^2}\)는 표본 데이터의 분산이다.

\[ \hat\mu = \dfrac{1}{N}\sum_{i=1}^{N} x_i \]

\[ \hat{\sigma^2} = \dfrac{1}{N}\sum_{i=1}^{N} (x_i - \hat\mu)^2 \]

위 식에서 \(x_i\)는 \(N\)개의 표본 데이터 각각의 값을 나타낸다.

이 모수 추정치는 각각 numpy 패키지의 mean, var 함수로 계산할 수 있다. 예를 들어 붓꽃 꽃받침의 길이의 평균과 분산은 5.006, 0.121766다.

mu_hat, sigma2_hat = np.mean(setosa_sepal_length), np.var(setosa_sepal_length)
mu_hat, sigma2_hat

(5.006, 0.12176400000000002)

scipy 패키지의 stats 서브패키지에서 제공하는 norm 클래스의 fit 정적 메서드를 사용해서 모수 추정을 할 수도 있다. 이 방법을 사용하면 정규분포 뿐 아니라 뒤에서 설명할 임의의 수치형 확률분포함수의 모수를 추정할 수 있다는 장점이 있다. 다만 norm 클래스에서는 분산이 아닌 분산의 제곱근 값인 표준편차(standard deviation)를 계산한다.

mu_hat, sigma_hat = norm.fit(setosa_sepal_length)
mu_hat, sigma_hat

(5.006, 0.3489469873777391)

sigma_hat ** 2

0.12176400000000001

사실 seaborn 패키지의 distplot 함수는 이 방법으로 모수를 추정하여 추정된 확률분포함수의 모양을 히스토그램에 겹쳐 그리는 기능을 제공하고 있다. distplot 함수의 fit 인수에 norm 클래스를 대입하면 된다.

sns.distplot(setosa_sepal_length, fit=norm, kde=False)
plt.title("붓꽃 꽃받침의 길이 분포와 추정된 정규분포함수의 확률밀도함수")
plt.show()

2.6.1.6. 정규분포의 모수검정

이렇게 추정한 정규분포의 모수값이 특정한 값과 같은지 아닌지를 검정을 통해 검증할 수 있다. 정규분포의 평균 모수를 검정할 때는 단일표본 t검정(one sample t-test)이라는 방법을 사용한다. 이 이름은 뒤에서 설명할 스튜던트t분포에서 나온 이름이다.

예를 들어 setosa종 붓꽃 꽃받침의 평균값은 5.06이 나왔지만 이 값을 그냥 5라고 할 수 있는지 알아보려면 다음과 같은 귀무가설 \(H_0\)에 대한 t검정을 하면 된다.

\[ H_0 : \mu = 5 \]

scipy 패키지의 stats 서브패키지는 t검정을 위한 ttest_1samp 함수를 제공한다. 이 함수에 데이터 값과 비교모수값을 인수로 넣으면 된다.

from scipy.stats import ttest_1samp

ttest_1samp(setosa_sepal_length, 5)

Ttest_1sampResult(statistic=0.12036212238318053, pvalue=0.9046884777690936)

유의확률이 90%가 넘으므로 모수값을 5라고 해도 문제가 없다고 할 수 있다.

독립표본 t검정(independent two sample t-test)을 사용하면 두 개의 정규분포의 모수값이 같은지 아닌지 비교할 수도 있다. 이 때는 같은 패키지의 ttest_ind 함수를 사용한다.

예를 들어 setosa 종의 꽃받침 길이는 5.06이었지만 versicolor 종의 꽃받침 길이는 5.936이다.

vericolor_sepal_length = iris[iris.Species == "versicolor"]["Sepal.Length"].values

np.mean(vericolor_sepal_length)

5.936

하지만 시각적 분포도에서 볼 수 있듯이 이 두 종의 꽃받침 길이 분포에는 겹치는 부분이 많다. 따라서 원래 두 종의 꽃받침 길이의 모수는 같지만 표본 데이터의 오차로 인해 이렇게 두 값이 다른 결과가 나왔을 수도 있다.

sns.histplot(x="Sepal.Length", hue="Species", element="poly",
             data=iris[iris.Species.isin(["setosa", "versicolor"])])
plt.show()

이 두 종의 원래 모수는 같았지만 표본 데이터의 오차에 의해 이렇게 다른 값이 나온건지 아닌지를 알아내려면 다음 코드를 사용한다. ttest_ind 함수에서 equal_var 인수는 두 표본데이터의 분산 모수가 같을 때는 True, 아니면 False를 넣는다. 일단 여기에서는 두 분산 모수가 같은지 모르므로 False 값을 넣는다.

from scipy.stats import ttest_ind

ttest_ind(setosa_sepal_length, vericolor_sepal_length, equal_var=False)

Ttest_indResult(statistic=-10.52098626754911, pvalue=3.746742613983842e-17)

유의확률이 \(3.7\times 10^{-17}\)로 두 종의 꽃받침 길이는 다르다는 것을 확인할 수 있다.

앞에서는 두 표본 데이터의 분산 모수가 다르다고 가정햐고 검정을 실시했는데 다음 함수를 사용하면 두 표본 데이터의 분산 모수가 같은지 다른지에 대한 검정을 실시할 수 있다.

• 바틀렛 검정 : scipy.stat.bartlett

• 플리그너 검정 : scipy.stat.fligner

• 레빈 검정 : scipy.stat.levene

이 함수들을 이용하여 setosa 종과 versicolor 종의 꽃받침 길이의 분산이 같은지 다른지 알아보자

from scipy.stats import bartlett, fligner, levene

bartlett(setosa_sepal_length, vericolor_sepal_length)

BartlettResult(statistic=6.891726740802407, pvalue=0.008659557933880048)

fligner(setosa_sepal_length, vericolor_sepal_length)

FlignerResult(statistic=7.394812196534449, pvalue=0.006541225398265243)

levene(setosa_sepal_length, vericolor_sepal_length)

LeveneResult(statistic=8.172720533728683, pvalue=0.005195521631017526)

3가지 검정의 결과가 모두 유의확률 1% 이내로 나왔으므로 두 종의 꽃받침 길이의 분산은 다르다는 것을 알 수 있다.

2.6.2. 로그정규분포

어떤 데이터는 그 자체로는 정규분포가 아니지만 로그 함수를 적용하면 정규분포가 되는 경우가 있다. 이러한 분포를 로그정규분포(lognormal distribution)라고 한다. 금액, 수량, 가격 등 양수 값만 존재하고 0이나 음수값이 존재하지 않는 데이터에서 이런 경우가 종종 발생한다.

2.6.2.1. 로그정규분포의 확률분포함수

확률변수 \(x\)가 정규분포를 따르는 경우 다음과 같이 표시한다.

\[ \log(x) \sim \mathcal{N}(x) \]

위 식에서 \(\mathcal{N}(x)\)은 정규분포의 확률밀도함수를 가리키는 기호다.

2.6.2.2. 로그정규분포의 시뮬레이션

scipy 패키지의 stats 서브패키지는 로그정규분포를 시뮬레이션할 수 있는 lognorm 클래스를 제공한다. lognorm 클래스의 인수는 loc와 s가 있다. loc 인수로 모수 \(\mu\)를 설정하고 s 인수로 모수 \(\sigma^2\)의 제곱근값 \(\sigma\)를 설정한다.

from scipy.stats import lognorm

mu = 0
s = 1
x = lognorm(loc=mu, s=np.sqrt(sigma2))
x

<scipy.stats._distn_infrastructure.rv_frozen at 0x138757370>

pdf 메서드를 사용하면 확률밀도함수를 계산할 수 있다. 인수로는 확률변수가 가질 수 있는 값을 넣는다.

xx = np.linspace(1e-10, 5, 100)
plt.plot(xx, x.pdf(xx))
plt.title("정규분포의 확률밀도함수")
plt.show()

다른 확률분포와 마찬가지로 rvs 메서드를 사용하여 표본측정 시뮬레이션을 할 수 있다.

import numpy as np

np.random.seed(0)

samples = x.rvs(100)

시뮬레이션된 표본값의 히스토그램을 그리면 다음과 같다.

import seaborn as sns

sns.displot(samples)
plt.xlim(0, 15)
plt.show()

2.6.2.3. 로그정규분포의 예

팁 데이터의 total_bill 데이터는 로그정규분포다. 이 데이터는 지불금액이므로 0이나 음수값이 존재하지 않고 그 자체는 정규분포를 띄지 않는다.

total_bill = tips["total_bill"]

sns.distplot(total_bill)
plt.show()

정규성 검정을 실시하면 유의확율이 \(1.595 \times 10^{-10}\)으로 정규분포가 아니다.

omni_normtest(total_bill)

NormaltestResult(statistic=45.11781912347332, pvalue=1.5951078766352608e-10)

하지만 이 값에 로그를 취하면 다음과 같이 정규분포 모양이 된다.

log_total_bill = np.log(total_bill)

sns.distplot(log_total_bill, fit=norm, kde=False)
plt.show()

정규성 검정을 실시하면 유의확률이 26.77%로 정규분포다.

omni_normtest(log_total_bill)

NormaltestResult(statistic=2.635400431729126, pvalue=0.26775036247744166)

2.6.2.4. 로그정규분포의 모수추정

로그정규분포의 모수는 다음과 같이 추정한다.

s, loc, scale = lognorm.fit(total_bill)
s, loc, scale 

(0.39848779123634986, -1.5965087437786663, 19.74259954314551)

팁 데이터의 total_bill 데이터에 대해 추정된 모수값으로 확률밀도함수를 그려서 원래 데이터의 히스토그램과 비교하면 다음과 같다.

x = lognorm(s=s, loc=loc, scale=scale)
xx = np.linspace(1e-10, 50)
sns.histplot(total_bill, stat='density')
plt.plot(xx, x.pdf(xx))
plt.show()

2.6.3. 스튜던트t분포

스튜던트t분포는 정규분포와 유사하지만 양 끝단의 꼬리 부분 데이터가 정규분포보다 많은 분포다.

2.6.3.1. 스튜던트t분포의 확률분포함수

확률변수 \(x\)가 정규분포를 따르는 경우 다음과 같이 표시한다.

\[ x \sim t(x) \]

위 식에서 \(t(x)\)은 스튜던트t분포의 확률밀도함수를 가리키는 기호다.

스튜던트t분포의 확률밀도함수 \(t(x)\)의 수식은 다음과 같다.

\[ t(x;\mu, \lambda, \nu) = \frac{\sqrt{\lambda}}{\sqrt{\nu\pi}}\frac{\Gamma\left(\frac{\nu+1}{2}\right)} {\Gamma\left(\frac{\nu}{2}\right)} \left(1+\lambda\frac{(x-\mu)^2}{\nu} \right)^{-\frac{\nu+1}{2}} \]

이 식에서 \(\lambda\)는 정규분포의 정밀도 \(\sigma^{-2}\)에 대응하는 개념이고 \(\Gamma(x)\)는 감마(Gamma) 함수라는 특수 함수다.

\[ \Gamma(x) = \int_0^\infty u^{x-1} e^{-u} d \]

스튜던트t분포는 양끝쪽의 꼬리부분의 두께를 조정할 수 있는 자유도(degree of freedom) 모수 \(\nu\)가 있다. 이 값은 2 이상의 자연수를 사용한다. 모수 \(\nu\)가 1인 경우는 코시분포(Cauchy distribution)라고 한다. 코시분포에서 양수인 부분만 사용하는 경우에는 하프코시분포(Half-Cauchy distribution)라고 부른다.

2.6.3.2. 스튜던트t분포의 시뮬레이션

scipy 패키지의 stats 서브패키지는 로그정규분포를 시뮬레이션할 수 있는 t 클래스를 제공한다. t 클래스의 인수는 df와 loc, scale 3가지 인데 df는 자유도, loc는 평균, scale은 표준편차를 설정한다. 메서드의 사용 방법은 정규분포와 같다.

from scipy.stats import t

df = 2
mu = 0
sigma2 = 1
x = t(df=df, loc=mu, scale=np.sqrt(sigma2))

xx = np.linspace(-5, 5, 100)
plt.plot(xx, x.pdf(xx))
plt.title("스튜던트t분포 확률밀도함수")
plt.show()

자유도 모수의 값이 증가하면 스튜던트t분포의 모양이 점점 정규분포와 닮아간다. 일반적으로 자유도가 30 이상이면 정규분포라고 본다.

for df in [1, 2, 5, 10, 20]:
    x = t(df=df)
    plt.plot(xx, x.pdf(xx), label=("스튜던트t분포(dof=%d)" % df))
plt.plot(xx, norm().pdf(xx), label="정규분포", lw=5, alpha=0.5)
plt.title("자유도에 따른 스튜던트t분포의 변화")
plt.xlabel("표본값")
plt.ylabel("확률밀도함수")
plt.legend()
plt.show()

2.6.3.3. 스튜던트t분포의 예

다음은 삼성전자 주식의 20년간의 주가수익률 데이터다.

import pandas_datareader.data as web

prices = web.DataReader('005930', 'naver', start='2001-1-1', end='2021-12-31')["Close"].astype(int)
prices

Date
2001-01-02     3350
2001-01-03     3370
2001-01-04     3780
2001-01-05     3900
2001-01-08     3880
              ...  
2021-12-24    80500
2021-12-27    80200
2021-12-28    80300
2021-12-29    78800
2021-12-30    78300
Name: Close, Length: 5191, dtype: int64

이 주식의 주가 수익률은 다음과 같이 구할 수 있다.

returns = prices.pct_change().dropna()

sns.distplot(returns)
plt.show()

주가수익률 데이터의 분포는 얼핏보면 정규분포와 유사한 형태를 가진다. 하지만 QQ플롯을 그려보면 직선이 아니라 양측 끝에서 휘어지는 모습을 볼 수 있다.

sp.stats.probplot(returns, plot=plt)
plt.show()

정규성 검정을 해도 유의확률이 \(4\times 10^{-103}\)으로 분포가 정규분포가 아니라는 것을 확인할 수 있다.

omni_normtest(returns)

NormaltestResult(statistic=471.43954606556974, pvalue=4.2481837795163223e-103)

이 분포가 정규분포를 따르지 않는 이유는 양 끝단의 데이터가 정규분포보다 많기 때문이다. 분포함수에서 꼬리부분이 길거나 두껍다는 뜻에서 이 현상을 롱테일(long tail) 혹은 팻테일(fat tail) 현상이라고 부른다. 위에서 예시로 든 주가 수익률의 경우 만약 정규분포를 따른다면 30년 기간동안 일수익률 7%를 넘는 수익률은 거의 발생하지 않았어야 했다. 하지만 실제로는 22번이나 발생했다.

len(returns[returns > 0.07])

22

히스토그램을 정규분포와 비교하여 삺펴보아도 정규분포의 꼬리 바깥쪽에 아웃라이어 데이터가 많이 존재하는 것을 볼 수 있다.

sns.distplot(returns, fit=norm, kde=False, rug=True, rug_kws={"height": 0.1,})
plt.xlim(0, 0.15)
plt.show()

2.6.3.4. 스튜던트t분포의 모수추정

스튜던트t분포의 모수는 다음과 같이 추정한다.

df, loc, scale = t.fit(returns)
df, loc, scale 

(4.6818243426797395, 0.0004703968628922298, 0.01567026026775519)

주가 수익률 데이터에 대해 추정된 모수값으로 확률밀도함수를 그려서 원래 데이터의 히스토그램과 비교하면 다음과 같다.

x = t(df=df, loc=loc, scale=scale)
xx = np.linspace(-0.1, 0.1)
sns.histplot(returns, stat='density')
plt.plot(xx, x.pdf(xx))
plt.show()