A B C d E F G H I K L M N O P Q R S T U V W Z A statsmodels.ap.stats.anova_lm(x) statsmodels.formula.api.ols 에 의해 생성되는 모형 즉, 클래스 인스턴스(x)를 인수로 받아 anova를 실행합니다. np.argsort(x, axis=-1, kind=None) 객체 x를 정렬할 경우 각 값에 대응하는 인덱스를 반환합니다. Axis는 기준 축을 지정하기 위한 매개변수로서 정렬의 방향을 조정할 수 있음(-1은 기본값으로 마지막 축) pandas.Series.autocorr(lag=1) lag에 전달한 지연수에 따른 값들 사이의 자기상관을 계산 B scipy.stats.bernoulli(x, p) 베르누이분포에 관련된 통계량을 계산하기 위한 클래스를 생성합니다. x: 랜덤변수 p: 단일 시행에서의 확률 scipy.stats.binom(x, n, p) 이항분포에 관련된 통계량을 계산하기 위한 클래스를 생성합니다. x: 랜덤변수 n: 총 시행횟수 p: 단일 시행에서의 확률 C scipy.stats.chi2.pdf(x, df, loc=0, scale=1) 카이제곱분포의 확률밀도함수를 계산 $$f(x, k) =\frac{1}{2^{\frac{k}{2}−1}Γ(\frac{k}{2})}x^{k−1}\exp\left(−\frac{x^2}{2}\right)$$ x: 확률변수 df: 자유도 pd.concat(objs, axis=0, join=’outer’, …) 두 개이상의 객체를 결합한 새로운 객체를 반환. objs: Series, DataFrame 객체. Axis=0은 행단위 즉, 열 방향으로 결합, Axis=1은 열단위 즉, 행 방향으
내용
정규성검정(Normality Test)
중심극한정리에 의해 자료의 샘플 수가 증가할수록 정규분포에 근접합니다. 특히 표본평균들의 분포인 표본분포는 정규분포에 부합합니다. 그러나 평균이 아닌 원시 데이터(raw data)의 경우 정규분포에 부합여부가 중요한 경우가 있습니다. 예를 들어 회귀분석의 경우 관찰값과 회귀 모형에 의한 예측값들의 차이를 잔차(residual라고 하는데 잔차가 정규분포에 부합한다는 가정하에 실시됩니다. 그 가정에 부합하는가의 여부가 성립된 모형의 적합도를 결정할 판단근거가 됩니다.
정규성 검정은 다음 방법을 사용합니다.
- Quantile-Quantile plot: 시각적 분석에 의한 결정
- Shaprio-Wilks test: 표본수(n < 2000)인 경우에 주로 사용
- Kolmogoroves-Smrinov test: n > 2000인 경우에 사용
Q-Q plot
Q-Q(사분위수) plot은 두 자료들을 분위수로 구분한 후 동일한 분위수 값들에 대해 작성한 도표로서 두 그룹의 분포를 비교하는 방법으로 역누적분포에 의해 설명됩니다. 역누적분포(Inverse cumulative distribution)는 누적분포의 역함수입니다. 예를 들어 표준정규분포의 누적분포와 역누적분포는 그림 1과 같이 나타낼 수 있습니다.
import numpy as np import pandas as pd from scipy import stats import matplotlib.pyplot as plt from scipy import special from sympy import * import FinanceDataReader as fdr
rng=np.linspace(-3, 3, 1000) pdf=[stats.norm.pdf(i) for i in rng] cdf=[stats.norm.cdf(i) for i in rng] inverCdf=[stats.norm.ppf(i) for i in cdf]
plt.figure(figsize=(10, 4))
plt.subplot(1,2,1)
plt.title("(A) PDF & CDF", size=12, weight="bold")
plt.plot(rng, pdf, label="PDF")
plt.plot(rng, cdf, label="CDF")
plt.xlabel("X", size=12, weight="bold")
plt.ylabel("Probability", size=12, weight="bold")
plt.legend(loc="best")
plt.subplot(1,2,2)
plt.title("(B) Inverse CDF = I(X)", size=12, weight="bold")
plt.plot(cdf, inverCdf, label="Inverse CDF")
plt.xlabel("CDF", size=12, weight="bold")
plt.ylabel("Probability", size=12, weight="bold")
plt.legend(loc="best")
plt.show()
어떤 두 그룹 A, B가 동일한 분포를 가진다면 각각의 역누적분포의 p값은 다음의 선형관계가 성립될 것입니다.
$$\begin{align}&I_A(x) = mI_B(x)+n\\ & m,n \in \mathbb{R} \end{align}$$위 식을 만족하는 두 그룹 중 B가 정규분포를 따른다면 A는 정규분포를 따른다고 할 수 있습니다. 물론 이 관계는 정성적 분석방법으로 부합 정도를 수치로 나타낼 수 없지만 대략적인 정규성을 판단할 수 있습니다.
다음은 일정구간의 랜덤수에 대한 Q-Q plot을 작성한 것입니다. Q-Q plot은
scipy.stats 클래스의 probplot()함수에 의해 작성합니다. 이 함수는
분석 대상의 표준화된 값, 원래 값, 그래프, 회귀선의 기울기와 편차(절편), 그리고
결정계수($R^2$)을 나타냅니다. 결정계수는 상관계수를 제곱한 것으로 [0, 1]범위의
값을 가집니다. 다음 코드의 결과 그림에서 회귀선과 완전부합하는 경우 1, 어떤
관계도 파악할 수 없는 경우는 0입니다. 그러므로 이 결과에 의해 정량적 판단이
가능해집니다.
x=np.random.rand(1000)
val, para=stats.probplot(x, plot=plt)
print(f'표준화된 값: {np.around(val[0][:3], 3)}
원시데이터: {np.around(val[1][:3], 3)}')
표준화된 값: [-3.198 -2.932 -2.784] 원시데이터: [0.003 0.006 0.006]
print(f'회귀선의 기울기: {round(para[0],3)}, 절편: {round(para[1],3)}, 결정계수: {round(para[1],3)}')
회귀선의 기울기: 0.283, 절편: 0.501, 결정계수: 0.501
예 1)
일정기간의 ksopi지수와 원-달러 환율 자료의 일일 종가(Close)에 대해 Q-Q
plot에 의한 정규성 검정을 실시해봅니다.
st=pd.Timestamp(2020,1,3)
et=pd.Timestamp(2022, 1, 14)
ks=fdr.DataReader('KS11', st, et)["Close"]
usKr=fdr.DataReader("USD/KRW", st, et)["Close"]
#kospi 지수 plt.figure(dpi=100) val, para=stats.probplot(ks, plot=plt, rvalue=True) np.around(para, 3)
array([4.653650e+02, 2.672327e+03, 9.470000e-01])
#원-달러 plt.figure(dpi=100) val3, para3=stats.probplot(usKr, plot=plt, rvalue=True) np.around(para3,3)
array([4.039100e+01, 1.162405e+03, 9.850000e-01])
shapiro-Wilk test
표본 $x_1, x_2, \cdots, x_n$이 정규분포에 부합성 여부를 검정하는 것으로 식 1의 W통계량에 대해 귀무가설(H0: 정규분포를 따릅니다.)에 대한 검정을 실시합니다.
$$\begin{align}\tag{1}&W=\frac{\left(\sum^n_{i=1} a_ix_{(i)}\right)^2}{\sum^n_{i=1} (x_i-\bar{x})^2}\\ &x_{(i)}=x_{n+1-i}-x_i\\ &n: \text{샘플수} \end{align}$$위 식에서 $a_i$는 두 값의 차이에 대한 가중치입니다. 그 가중치들은 정렬된 자료의 평균, 표준편차 등 통계량을 기준으로 산출된 상수로서 shapiro-Wilk table에서 결정할 수 있습니다. 또한 m은 샘플수를 기준으로 산출되는 값으로 다음과 같습니다.
scipy라이브러리의 stats.shapiro() 함수에 의해 검정할 수 있습니다.
이 함수는 W값을 반환하는 것으로 shapiro-Wilk table을 대신할 수 있습니다.
위 W 통계량 계산에서 $x_{(i)}$에서 i는 [0, m] 사이의 정수값입니다. 즉 W값은 전체 퍼짐의 정도에서 각각의 작은값과 큰값의 차이의 비를 나타낸 값입니다. Shapiro-Wilk 검정의 W는 다음 과정으로 계산합니다.
- data 정렬
- SS 계산 $$SS=\sum^n_{i=1} (x_i-\bar{x})^2$$
- W의 분자인 b를 계산 $$b=\sum^m_{i=1}a_i(x_{n+1-i}-x_i)$$
- 검정 통계량(W) 계산 $$W=\frac{b^2}{SS}$$
- shapiro-Wilk table를 기준으로 p-value를 산출. shapiro-Wilk table은 샘플수에 대한 가중치($a_i$)와 특정 유의확률에 대응하는 W값을 계산한 표입니다.
예 2)
다음 나이에 대한 자료가 정규분포에 부합 여부에 대해 shapiro-Wilk 검정을
실시해 봅니다.
age=np.array([65,61,63,86,70,55,74,35,72,68,45,58]) ageSort=np.sort(age) da=pd.DataFrame([age, ageSort], index=['age', 'ageSort']) da
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | |
| age | 65 | 61 | 63 | 86 | 70 | 55 | 74 | 35 | 72 | 68 | 45 | 58 |
| ageSort | 35 | 45 | 55 | 58 | 61 | 63 | 65 | 68 | 70 | 72 | 74 | 86 |
위 자료는 샘플수(n) 12이므로 $m=\frac{12}{2}=6$로서 6개의 가중치가 존재합니다. 각 가중치를 shapiro-Wilk table에서 결정하면 다음과 같습니다.
n=len(age) m=np.array([n/2 if n%2==0 else (n-1)/2], dtype=np.int16) m
array([6], dtype=int16)
ss=np.sum((ageSort-np.mean(ageSort))**2) round(ss, 3)
2008.667
a=np.array([0.5475,0.3325,0.2347,0.1586,0.0922,0.0303])#shapiro-Wilk table에서 발견할 수 있음 age1=[] for i in range(m[0]): age1.append(ageSort[n-(i+1)]-ageSort[i]) age1
[51, 29, 17, 12, 7, 2]
b=np.sum(a*age1) b
44.1641
w=b**2/ss round(w, 3)
0.971
shapiro-Wilk table에서 n=12인 경우 유의확률과 W값은 다음과 같습니다. 결정된
값을 기준으로 선형적으로 산출된 w에 대한 유의확률을 추정하면 다음과 같습니다.
우선 다음 값들에서 직선식을 결정합니다. 이 직선식은 수기로 충분히 작성할 수
있습니다. 다음은 sympy 모듈의 solve()함수를 사용한 것입니다.
tab=pd.DataFrame([[0.5, 0.943],[0.9, 0.973]], columns=['유의확률','W']) tab
| 유의확률 | W | |
|---|---|---|
| 0 | 0.5 | 0.943 |
| 1 | 0.9 | 0.973 |
tab1=tab.values slope=(tab1[0,0]-tab1[1,0])/(tab1[0,1]-tab1[1,1]) round(slope,3)
13.333
from sympy import*
x, bias=symbols('x bias')
eq=slope*x+bias
eq
bias+13.3333333333333x
eq1=eq.subs(x,tab1[0,1]) eq1
bias+12.5733333333333
bias1=solve(Eq(eq1, tab1[0,0]), bias) bias1[0].evalf(4)
−12.07
linear=eq.subs(bias, bias1[0]) N(linear,4)
13.33𝑥−12.07
위에서 산출된 w값을 직선식에 적용하면 대략적인 유의확률을 반환합니다.
p=linear.subs(x, w) round(p,3)
0.874
위 계산과정은 scpy.stats.shapiro()함수로 대신할 수 있습니다. 이
함수는 w값과 p-value를 반환합니다.
statics, pVal=stats.shapiro(age)
print(f'통계량: {round(statics, 3)}, p-value: {round(pVal, 3)}')
통계량: 0.971, p-value: 0.922
위 결과의 p 값은 이전에 외삽(interpolation)으로 계산한 대략적인 p값과 차이를 보입니다. 그러나 두 값 모두 유의수준 0.05 또는 0.1보다 크므로 위 자료 age는 정규분포를 따른다는 귀무가설을 기각할 수 없습니다. 이 결과를 QQ plot 분석과 비교해 보면 다음과 같습니다.
val, para=stats.probplot(age, plot=plt, rvalue=True)
print(f'회귀계수: {np.around(para[0], 3)}\n편차: {np.around(para[1], 3)}\n결정계수($R^2$): {np.around(para[2], 3)}')
회귀계수: 14.264 편차: 62.667 결정계수($R^2$): 0.979
예 3)
예 2에서 사용된 kospi에 대해 shapiro-Wilk 검정을 실시해 봅니다.
statics, pVal=stats.shapiro(ks.values)
print(f'통계량: {round(statics, 3)}, p-value: {round(pVal, 3)}')
통계량: 0.895, p-value: 0.0
유의 수준 0.05에서 위 결과는 Q-Q plot과는 달리 정규분포에 부합하는다는 귀무가설을 채택할 수 없습니다.
Kolmogorov-Smirnov Test
Kolmogorov-Smirnov 검정 (K-S test)은 표본이 특정 분포를 가진 모집단에서 추출되었는지 결정하는 데 사용됩니다. K-S 검정은 경험적 분포 함수 (ECDF)를 기반으로 하며 N 개의 정렬된 데이터 포인트 $Y_1, Y_2, \cdots, Y_N$이 주어지면 ECDF는 식 2와 같이 정의됩니다.
$$\begin{equation}\tag{2}\text{ECDF}=\frac{n(i)}{N}\end{equation}$$여기서 n(i)는 $Y_i$보다 적은 포인트 수이고 $Y_i$는 오름차순으로 정렬됩니다. 정렬된 각 데이터 포인트의 값에서 $\displaystyle \frac{1}{N}$ 씩 증가하는 단계 함수입니다.<\p>
그림 2는 100 개의 정규 난수에 대해 정규 누적 분포 함수와 경험적 분포 함수의 플롯입니다. K-S 테스트는 이 두 곡선 사이의 최대 거리를 기반으로합니다.
np.random.seed(3) da=np.sort(np.random.randn(100)) np.around(da[:3], 3)
array([-2.916, -2.419, -2.248])
ECDF=[]
N=len(da)
for i in range(N):
ECDF.append((i+1)/N)
ECDF[:3]
[0.01, 0.02, 0.03]
NCDF=[stats.norm.cdf(i) for i in da] np.around(NCDF[:3], 3)
array([0.002, 0.008, 0.012])
plt.figure(dpi=100)
plt.plot(da, ECDF, color="blue", label="ECDF")
plt.plot(da, NCDF, color="red", label="NCDF")
plt.xlabel("X", size="13")
plt.ylabel("Cumul. Prob.", size="13")
plt.legend(loc="best")
plt.title("100 random number", size="14")
plt.show()
이 검정은 몇 가지 중요한 제한 사항이 있습니다.
- 연속 분포에만 적용
- 꼬리보다 분포의 중심 근처에서 더 민감
- 분포를 완전히 지정해야한다는 것입니다. 즉, 데이터에서 위치, 축척 및 모양 매개 변수를 추정하면 K-S 테스트의 임계 영역이 더 이상 유효하지 않습니다. 일반적으로 시뮬레이션에 의해 결정되어야합니다.
K-S 검정의 가설은 다음과 같습니다.
- H0: 데이터는 특정한 분포를 따릅니다.
- H1: 데이터는 특정한 분포를 따르지 않습니다.
F는 연속변수의 이론적 누적 분포 함수로서 평균(location), 편차정도(scale), 형상(shape) 매개변수는 데이터로부터 추정된 것이 아닙니다. 유의수준을 기준으로 통계량의 기각 여부를 결정하기 위해서 Kolmogorov table의 임계값(critical value)를 적용합니다. 임계값은 K-S table에 의해 샘플의 규모와 유의수준에 의해 결정되지만 샘플의 크기가 50 이상일 경우 다음과 같이 계산됩니다.
| n|$\alpha$ | 0.001 | 0.01 | 0.02 | 0.05 | 0.1 | 0.15 | 0.2 |
|---|---|---|---|---|---|---|---|
| n>50 | $\frac{1.95}{n}$ | $\frac{1.63}{n}$ | $\frac{1.52}{n}$ | $\frac{1.36}{n}$ | $\frac{1.22}{n}$ | $\frac{1.14}{n}$ | $\frac{1.07}{n}$ |
검정통계량 D와 임계값의 비교에 의해 다음과 같은 결론을 유도됩니다.
$$D > \text{임계값} \rightarrow H0: \text{기각}$$
예 4)
예 2에서 사용한 kospi 지수의 일일변동률 자료에 대해 K-S 검정 합니다.
다음 과정으로 분석합니다.
- 자료의 표준화
- 오름차순으로 정렬
- 표준정규 누적분포와 경험적 누적분포를 계산
- 검정 통계량 계산
- 임계값 계산
emCDF=[] # 경험적누적분포
N=len(da) # 자료 크기
for i in range(N):
emCDF.append((i+1)/N)
np.around(emCDF[:3],3)
array([0.01, 0.02, 0.03])
norCDF=[stats.norm.cdf(i) for i in da] # 표준정규 누적분포f np.around(norCDF[:3],3)
array([0.002, 0.008, 0.012])
plt.figure(dpi=100)
plt.plot(da, emCDF, color="blue", label="emCDF")
plt.plot(da, norCDF, color="red", label="norCDF")
plt.xlabel("X", size="13")
plt.ylabel("Cumul. Prob.", size="13")
plt.legend(loc="best")
plt.title("Kospi K-S Test", size="14")
plt.show()
D=np.max(abs(np.array(emCDF)-np.array(norCDF))) # K-S 통계량
print(f'통계량: {round(D,3)}')
통계량: 0.097
#critical values(임계값) Kolmogorov table로 부터 정해진 값
crit=1.36/np.sqrt(N)
print(f'임계값: {round(crit,3)}')
임계값: 0.136
위의 검정통계량(D)는 임계값에 비해 작습니다. 그러므로 귀무가설을 기각할 수 없습니다. 이 검정은 scipy.stats.kstest()함수를 사용합니다. 결과로 반환된 p-value는 유의수준 0.05 또는 0.01 크므로 귀무가설을 기각할 수 없습니다.
statics, pVal=stats.kstest(da, 'norm')
print(f'통계량: {round(statics, 3)}, p-value:{round(pVal, 3)}')
통계량: 0.097, p-value:0.282
위 결과는 shapiro 검정의 결과와는 다르지만 Q-Q plot와는 동일한 결론을 보입니다. shapiro 검정의 경우는 다른 검정방법에 비해 매우 보수적이며 작은 규모의 자료에서 유용하지만 일정 규모 이상의 경우 크기에 민감합니다.
Jarque-Bera test
통계에서 Jarque-Bera 검정은 표본 데이터에 정규 분포와 일치하는 왜도 및 첨도가 있는지 여부에 대한 적합도(goodness-of-fit) 검정으로 검정통계량은 식 4와 같이 정의됩니다. 검정 통계량은 항상 음이 아닙니다. 0에서 멀다면 데이터가 정규 분포를 갖지 않는다는 신호입니다. 테스트 통계 JB는 다음과 같이 정의됩니다. $$\begin{align}&\tag{4}\text{JB}=\frac{n}{6}\left(\text{Skew}^2+\frac{1}{4}(\text{Kurt}-3)^2 \right)\\&n: \text{size of data}\\&\text{skew: 왜도}\\&\text{kurt: 첨도}\end{align}$$
왜도와 첨도는 3차 4차 중심모멘트로 계산됩니다.
데이터가 정규 분포에서 나온 경우 JB 통계량은 자유도가 2인 카이 제곱 분포를 점근적으로 가지므로 데이터가 정규 분포에서 왔다는 가설을 테스트하는 데 통계량을 사용할 수 있습니다. 정규분포의 왜도와 첨도는 각각 0과 3입니다. 그러므로 귀무가설은 다음과 같습니다.
식 4인 JB의 정의에서 알 수 있듯이 정규분포에 벗어나면 JB 통계량은 증가합니다. 작은 표본의 경우 카이 제곱 근사는 지나치게 민감하여 귀무 가설이 참일 때 종종 기각됩니다. 또한 p-값의 분포는 균일 분포에서 벗어나 특히 작은 p-값의 경우 오른쪽으로 치우친 단봉 분포가 됩니다. 이는 큰 제1종 오류율로 이어집니다.
JB 검정은 scipy.stats 모듈의 jarque_bera() 함수를 적용할 수 있습니다. 이 함수는 검정통계량과 p-value를 반환합니다.
예 5)
정규분포를 따르는 랜덤수에 대한 JB 검정을 실시해봅니다.
x=stats.norm.rvs(size=1000, random_state=0)
x1=np.sort(x)
y=[stats.norm.pdf(i) for i in x1]
plt.figure(dpi=100)
plt.hist(x, bins=30, density=True, alpha=0.3, label="hist")
plt.plot(x1, y, label="N(0, 1)")
plt.xlabel("data", size=12, weight="bold")
plt.ylabel("pdf, density", size=12, weight="bold")
plt.legend(loc="best", prop={'weight':"bold"})
plt.show()
jbtest=stats.jarque_bera(x) jbtest
Jarque_beraResult(statistic=0.28220016508625245, pvalue=0.8684023954281483)
위 결과에 의하면 유의수준 0.05에 비해 매우 높은 p값을 나타냅니다. 즉, 위 자료는 정규분포에 부합한다는 (JB ≈ 0) 귀무가설을 기각할 수 없습니다.
예 6)
랜덤수에 대한 JB 검정을 실시해봅니다.
x=np.random.rand(1000)
x1=np.sort(x)
y=[stats.norm.pdf(i) for i in x1]
plt.figure(dpi=100)
plt.hist(x, bins=30, density=True, alpha=0.3, label="hist")
plt.plot(x1, y, label="N(0, 1)")
plt.xlabel("data", size=12, weight="bold")
plt.ylabel("pdf, density", size=12, weight="bold")
plt.legend(loc="best", prop={'weight':"bold"})
plt.show()
stats.jarque_bera(x)
Jarque_beraResult(statistic=59.798510944727006, pvalue=1.0347278589506459e-13)
댓글
댓글 쓰기