均值向量的推断

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]]