import numpy as np
from operator import *
from cytoolz import *
from scipy import stats
def Hotelling(X, u0_, alpha):
n,p = X.shape
x_ = np.array(list(map(lambda x:reduce(add, X[:,x])/n, range(X.shape[1]))))
S = np.eye(p)
for i in range(p):
for j in range(p):
if i == j:
S[i, j] = reduce(add, map(lambda x: pow(x-x_[i], 2), X[:,i]))/(n-1)
else:
S[i, j] = reduce(add, map(lambda x: (x[0]-x_[i])*(x[1]-x_[j]), zip(X[:,i], X[:,j])))/(n-1)
S_r = np.linalg.inv(S)
T = n*(x_ - u0).dot(S_r).dot(x_.T - u0.T)
F = stats.f.isf(alpha, p, n-p)
print("u:")
print(x_)
print("S:")
print(S)
if T > F:
print(f"拒绝H0假设 T:{T:0.3f} > F:{F:0.3f}")
print('这个检验假设没通过啊')
else:
print(f"接受均值假设T:{T:0.3f} < F:{F:0.3f}")
X = np.array([
[3.7, 5.7, 3.8, 3.2, 3.1, 4.6, 2.4, 7.2, 6.7, 5.4, 3.9, 4.5, 3.5, 4.5, 1.5, 8.5, 4.5, 6.5, 4.1, 5.5],
[48.5, 65.1, 47.2, 53.2, 55.5, 36.1, 24.8, 33.1, 47.4, 54.1, 36.9, 58.8, 27.8, 40.2, 13.5, 56.4, 71.6, 52.8, 44.1, 40.9],
[9.3, 8.0, 10.9, 12.0, 9.7, 7.9, 14.0, 7.6, 8.5, 11.3, 12.7, 12.3, 9.8, 8.4, 10.1, 7.1, 8.2, 10.9, 11.2, 9.4]]).T
u0 = np.array([4, 50, 10])
Hotelling(X, u0, 0.1)
u:
[ 4.64 45.4 9.965]
S:
[[ 2.87936842 10.01 -1.80905263]
[ 10.01 199.78842105 -5.64 ]
[ -1.80905263 -5.64 3.62765789]]
拒绝H0假设 T:9.739 > F:2.437
这个检验假设没通过啊
def confidence_ellipse(X, alpha):
n,p = X.shape
x_ = np.array(list(map(lambda x:reduce(add, X[:,x])/n, range(X.shape[1])))) #椭球中心
S = np.eye(p)
for i in range(p):
for j in range(p):
if i == j:
S[i, j] = reduce(add, map(lambda x: pow(x-x_[i], 2), X[:,i]))/(n-1)
else:
S[i, j] = reduce(add, map(lambda x: (x[0]-x_[i])*(x[1]-x_[j]), zip(X[:,i], X[:,j])))/(n-1)
S_r = np.linalg.inv(S)
l, e = np.linalg.eig(S)
T = n*(x_ - u0).dot(S_r).dot(x_.T - u0.T)
F = stats.f.isf(alpha, p, n-p)
upper = list(map(lambda x:np.sqrt(x)*np.sqrt(p*(n-1)*F/(n*(n-p))), l))# 长轴长度
return x_, e, upper
confidence_ellipse(X, 0.1)
(array([ 4.64 , 45.4 , 9.965]),
array([[-0.05084144, -0.81748351, -0.57370364],
[-0.99828352, 0.02487655, 0.05302042],
[ 0.02907156, -0.57541452, 0.81734508]]),
[9.05067409589578, 0.7292367091847594, 1.3607856615453258])
def Joint_confidence_interval(X, alpha):
n,p = X.shape
x_ = np.array(list(map(lambda x:reduce(add, X[:,x])/n, range(X.shape[1])))) #椭球中心
S = np.eye(p)
for i in range(p):
for j in range(p):
if i == j:
S[i, j] = reduce(add, map(lambda x: pow(x-x_[i], 2), X[:,i]))/(n-1)
else:
S[i, j] = reduce(add, map(lambda x: (x[0]-x_[i])*(x[1]-x_[j]), zip(X[:,i], X[:,j])))/(n-1)
F = stats.f.isf(alpha, p, n-p)
off = list(map(lambda x:np.sqrt(p*(n-1)*F/(n-p))*np.sqrt(S[x,x]/n), range(p)))# 长轴长度
ran = list(map(lambda x:[x_[x]-off[x],x_[x]+off[x]], range(p)))
return ran
Joint_confidence_interval(X, 0.1)
def Single_confidence_interval(X, alpha):
n,p = X.shape
x_ = np.array(list(map(lambda x:reduce(add, X[:,x])/n, range(X.shape[1])))) #椭球中心
S = np.eye(p)
for i in range(p):
for j in range(p):
if i == j:
S[i, j] = reduce(add, map(lambda x: pow(x-x_[i], 2), X[:,i]))/(n-1)
else:
S[i, j] = reduce(add, map(lambda x: (x[0]-x_[i])*(x[1]-x_[j]), zip(X[:,i], X[:,j])))/(n-1)
t = stats.t.isf(alpha/2, n-1)
off = list(map(lambda x:t*np.sqrt(S[x,x]/n), range(p)))# 长轴长度
ran = list(map(lambda x:[x_[x]-off[x],x_[x]+off[x]], range(p)))
return ran
Single_confidence_interval(X ,0.1)
[[3.9839121839272464, 5.2960878160727525],
[39.93489498718458, 50.86510501281542],
[9.228578403535657, 10.701421596464343]]
def Bonferroni(X, alpha):
n,p = X.shape
x_ = np.array(list(map(lambda x:reduce(add, X[:,x])/n, range(X.shape[1])))) #椭球中心
S = np.eye(p)
for i in range(p):
for j in range(p):
if i == j:
S[i, j] = reduce(add, map(lambda x: pow(x-x_[i], 2), X[:,i]))/(n-1)
else:
S[i, j] = reduce(add, map(lambda x: (x[0]-x_[i])*(x[1]-x_[j]), zip(X[:,i], X[:,j])))/(n-1)
t = stats.t.isf(alpha/(2*p), n-1)
off = list(map(lambda x:t*np.sqrt(S[x,x]/n), range(p)))# 长轴长度
ran = list(map(lambda x:[x_[x]-off[x],x_[x]+off[x]], range(p)))
return ran
Bonferroni(X, 0.1)
[[3.7694351029105735, 5.510564897089425],
[38.14833571587469, 52.65166428412531],
[8.98783996801227, 10.94216003198773]]
