image.png

image.png

import numpy as np
from scipy import stats
x=np.array([[72,66,76,77],[60,53,66,63],[56,57,64,58],[41,29,36,38],
           [32,32,35,36],[30,35,34,26],[39,39,31,27],[42,43,31,25],[37,40,31,25],
                        [33,29,27,36],[32,30,34,28],[63,45,74,63],[54,46,60,52],
                        [47,51,52,43],[91,79,100,75],[56,68,47,50],[79,65,70,61],
                        [81,80,68,58],[78,55,67,60],[46,38,37,38],[39,35,34,37],
                        [32,30,30,32],[60,50,67,54],[35,37,48,39],[39,36,39,31],
                        [50,34,37,40],[43,37,39,50],[48,54,57,43]])
z=x.shape#(28,4)
n=z[0]
p=z[1]
dfchi=p*(p+1)*(p+2)/6#
q=np.eye(n)-1/n*np.ones((n,n))
s=1/n*np.dot(np.dot(x.T,q),x)
s_inv=np.linalg.inv(s)#矩阵求逆
g_matrix=np.dot(np.dot(np.dot(np.dot(q,x),s_inv),x.T),q)#矩阵相乘
print("based on skewness")
beta1hat=(sum(sum(g_matrix*g_matrix*g_matrix)))/(n*n) #矩阵对应元素相乘
print(beta1hat)
kappa1=n*beta1hat/6
print(kappa1)
pvalshew=1 - stats.chi2.cdf(kappa1,dfchi)
    #1-cdf('chi2',kappa1,cdfchi)
print(pvalshew)

print("based on kurtosis:")
beta2hat=(g_matrix*g_matrix).trace()/n
print(beta2hat)
kappa2=(beta2hat-p*(p+2))/np.sqrt(8*p*(p+2)/n)
print(kappa2)
pvalkurt = 2*(1 - stats.norm.cdf(abs(kappa2)))
    #2*(1-cdf('norm',abs(kappa2),0,1))
print(pvalkurt)

image.png