기본 콘텐츠로 건너뛰기

[matplotlib] 등고선(Contour)

[data analysis]shapiro-Wilk test

shapiro-Wilk test

관련된 내용

표본 x1, x2, …, xn이 정규분포에 부합성 여부를 검정하기 위해 식 1과 같이 계산되는 shapiro-Wilk통계량인 W를 사용하여 귀무가설(H0: 정규분포를 따릅니다.) 검정을 실시합니다.

W=(i=1maix(i))2i=1n(xix¯)2x(i)=xn+1ixin:표본의 크기m={n2forn=짝수n12forn=홀수 (식 1)

식 1에서 ai는 두 값 차이에 대한 가중치입니다. 그 가중치들은 정렬된 자료의 평균, 표준편차 등 통계량을 기준으로 산출된 상수로서 shapiro-Wilk table에서 결정할 수 있습니다.

식 1의 x(i)를 나타내는 i는 [0, m] 사이의 정수값입니다. 즉 W값은 전체 퍼짐의 정도에서 각각의 작은값과 큰값의 차이의 비를 나타낸 값입니다. Shapiro-Wilk 검정의 W는 다음 과정으로 계산합니다.

  1. data 정렬
  2. SS 계산
    • SS=i=1n(xix¯)2
  3. W의 분자인 b를 계산
    • b=i=1mai(xn+1ixi)
  4. 검정 통계량(W)계산
    • W=b2SS
  5. shapiro-Wilk table를 기준으로 p-value를 산출
    • shapiro-Wilk table은 샘플수에 대한 가중치(ai)와 특정 유의확률에 대응하는 W값을 계산한 표입니다.

이 검정은 scipy라이브러리의 stats.shapiro() 함수를 사용하여 실시할 수 있습니다.

예 1)

다음 나이에 대한 자료가 정규분포에 부합 여부에 대해 shapiro-Wilk 검정을 적용합니다.

age 65 61 63 86 70 55 74 35 72 68 45 58

위 자료를 내림차순으로 정리합니다.

age=np.array([65,61,63,86,70,55,74,35,72,68,45,58])
ageSort=np.sort(age)
print(ageSort)
[35 45 55 58 61 63 65 68 70 72 74 86]

위 자료는 샘플수(n)이 12이므로 m=122=6으로서 6개의 가중치가 존재합니다. 각 가중치를 shapiro-Wilk table에서 결정하고 통계량 W를 계산합니다.

n=len(age)
m=np.array([n/2 if n%2==0 else (n-1)/2], dtype=np.int16)
print(m)
[6]
ss=np.sum((ageSort-np.mean(ageSort))**2)
print("ss: %.3f" %ss)
ss: 2008.667

다음 코드의 객체 a는 shapiro-Wilk table로부터 결정되는 값들입니다.

a=np.array([0.5475,0.3325,0.2347,0.1586,0.0922,0.0303])
age1=[ageSort[n-(i+1)]-ageSort[i] for i in range(m[0])];age1
[51, 29, 17, 12, 7, 2]
b=np.sum(a*age1); b
44.1641
w=b**2/ss; w
0.9710260847046132

위 결과의 통계량에 대한 유의확률은 shapiro-Wilk table로부터의 값으로 외삽법을 사용하여 계산할 수 있습니다. 이 부분은 생략하며 stats.shapiro() 함수를 사용하여 유의확률을 확인합니다.

re=stats.shapiro(age)
print("통계량(w): %.3f, p-value: %.3f" %(re[0], re[1]))
통계량(w): 0.971, p-value: 0.922

유의수준 0.05 또는 0.1보다 크므로 위 자료 age는 정규분포를 따른다는 귀무가설을 기각할 수 없습니다. 이 결과에 대한 Q-Q plot은 그림 1과 같습니다.

그림 1. 자료 age에 대한 Q-Q plot.
plt.figure(figsize=(4,2))
stats.probplot(age, plot=plt, rvalue=True)
plt.show()

예 2)

다음은 일정한 기간의 kospi 지수와 kosdaq 지수의 일일 종가 자료입니다.

kospi kosdaq
0 2669.8 878.9
1 2607.3 871.6
2 2587.0 866.2
3 2578.1 878.3
4 2567.8 879.3

각 자료의 정규성을 파악하기 위해 shapiro-Wilk 검정을 실시해 봅니다.

위 자료를 생성하기 위한 코드입니다. 자료간의 스케일의 차이를 보정하기 위해 각 자료를 표준화하였습니다.

st=pd.Timestamp(2024,1,1)
et=pd.Timestamp(2024, 5, 30)
kos=fdr.DataReader("KS11",st, et)["Close"]
kq=fdr.DataReader("KQ11", st, et)["Close"]
kos1=(kos-kos.mean())/kos.std()
kq1=(kq-kq.mean())/kq.std()

위 코드로부터 작성된 자료 kos1과 ex1에 대한 Shapiro-Wilk 분석을 실행합니다.

kos_sh=stats.shapiro(kos1)
print("통계량: %.3f, p-value: %.3f" %(kos_sh[0], kos_sh[1]))
통계량: 0.921, p-value: 0.000
kq_sh=stats.shapiro(kq1)
print("통계량: %.3f, p-value: %.3f" %(kq_sh[0], kq_sh[1]))
통계량: 0.984, p-value: 0.250

유의수준 0.05를 기준으로 kospi는 매우 낮은 p-value를 보이므로 정규성을 따른다는 귀무가설을 기각할 수 있습니다. 반면에 kosdaq의 경우는 귀무가설을 기각할 수 없습니다.이 두 결과는 그림 2와 같이 시각화 할 수 있습니다.

그림 2. kospi와 kosdaq의 Q-Q plot.
plt.figure(figsize=(8, 3))
ax1 = plt.subplot(121)
res = stats.probplot(kos1, plot=plt, rvalue=True)
ax1.set_title("Q-Q plot: KOSPI")
ax2 = plt.subplot(122)
res = stats.probplot(kq1, plot=plt, rvalue=True)
ax2.set_title("Q-Q plot: kq")
ax2.set_ylabel("")
plt.show()

댓글

이 블로그의 인기 게시물

[Linear Algebra] 유사변환(Similarity transformation)

유사변환(Similarity transformation) n×n 차원의 정방 행렬 A, B 그리고 가역 행렬 P 사이에 식 1의 관계가 성립하면 행렬 A와 B는 유사행렬(similarity matrix)이 되며 행렬 A를 가역행렬 P와 B로 분해하는 것을 유사 변환(similarity transformation) 이라고 합니다. (1)A=PBP1P1AP=B 식 2는 식 1의 양변에 B의 고유값을 고려한 것입니다. (식 2)BλI=P1APλP1P=P1(APλP)=P1(AλI)P 식 2의 행렬식은 식 3과 같이 정리됩니다. det(BλI)=det(P1(APλP))=det(P1)det((AλI))det(P)=det(P1)det(P)det((AλI))=det(AλI)det(P1)det(P)=det(P1P)=det(I) 유사행렬의 특성 유사행렬인 두 정방행렬 A와 B는 'A ~ B' 와 같...

[sympy] Sympy객체의 표현을 위한 함수들

Sympy객체의 표현을 위한 함수들 General simplify(x): 식 x(sympy 객체)를 간단히 정리 합니다. import numpy as np from sympy import * x=symbols("x") a=sin(x)**2+cos(x)**2 a sin2(x)+cos2(x) simplify(a) 1 simplify(b) x3+x2x1x2+2x+1 simplify(b) x - 1 c=gamma(x)/gamma(x-2) c Γ(x)Γ(x2) simplify(c) (x2)(x1) 위의 예들 중 객체 c의 감마함수(gamma(x))는 확률분포 등 여러 부분에서 사용되는 표현식으로 다음과 같이 정의 됩니다. 감마함수는 음이 아닌 정수를 제외한 모든 수에서 정의됩니다. 식 1과 같이 자연수에서 감마함수는 factorial(!), 부동소수(양의 실수)인 경우 적분을 적용하여 계산합니다. (식 1)Γ(n)={(n1)!n:자연수0xn1exdxn:부동소수 x=symbols('x') gamma(x).subs(x,4) 6 factorial 계산은 math.factorial() 함수를 사용할 수 있습니다. import math math.factorial(3) 6 a=gamma(x).subs(x,4.5) a.evalf(3) 11.6 simpilfy() 함수의 알고리즘은 식에서 공통사항을 찾아 정리하...

sympy.solvers로 방정식해 구하기

sympy.solvers로 방정식해 구하기 대수 방정식을 해를 계산하기 위해 다음 함수를 사용합니다. sympy.solvers.solve(f, *symbols, **flags) f=0, 즉 동차방정식에 대해 지정한 변수의 해를 계산 f : 식 또는 함수 symbols: 식의 해를 계산하기 위한 변수, 변수가 하나인 경우는 생략가능(자동으로 인식) flags: 계산 또는 결과의 방식을 지정하기 위한 인수들 dict=True: {x:3, y:1}같이 사전형식, 기본값 = False set=True :{(x,3),(y,1)}같이 집합형식, 기본값 = False ratioal=True : 실수를 유리수로 반환, 기본값 = False positive=True: 해들 중에 양수만을 반환, 기본값 = False 예 x2=1의 해를 결정합니다. solve() 함수에 적용하기 위해서는 다음과 같이 식의 한쪽이 0이 되는 형태인 동차식으로 구성되어야 합니다. x21=0 import numpy as np from sympy import * x = symbols('x') solve(x**2-1, x) [-1, 1] 위 식은 계산 과정은 다음과 같습니다. x21=0(x+1)(x1)=0x=1or1x4=1의 해를 결정합니다. solve() 함수의 인수 set=True를 지정하였으므로 결과는 집합(set)형으로 반환됩니다. eq=x**4-1 solve(eq, set=True) ([x], {(-1,), (-I,), (1,), (I,)}) 위의 경우 I는 복소수입니다.즉 위 결과의 과정은 다음과 같습니다. x41=(x2+1)(x+1)(x1)=0x=±1,±1=±i,±1 실수...