Task4：数理统计

# 引入相关工具库
import numpy as np 
import pandas as pd 
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use("ggplot")
import warnings 
warnings.filterwarnings("ignore")
plt.rcParams['font.sans-serif']=['SimHei','Songti SC','STFangsong']
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
import seaborn as sns

频数直方图

x_samples = np.random.randn(1000)
plt.hist(x_samples, bins=10,color='blue',alpha=0.6)  # bins=10代表10根柱子
plt.xlabel("x")
plt.ylabel("频数 n")
plt.title("频数直方图")
plt.show()

频率直方图

x_samples = np.random.randn(1000)
plt.hist(x_samples, bins=10,color='blue',alpha=0.6,density=True)  # bins=10代表10根柱子
plt.xlabel("x")
plt.ylabel("频率 p")
plt.title("频率直方图")
plt.show()

从总体/总体的分布中抽取样本并计算样本均值和计算偏差

## （1）从总体中抽取样本
X = np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20])   # 假设总体为X
x_i = np.random.choice(X, 10, replace=False)  # 从总体X中抽取10个样本
x_mean = np.mean(x_i) # 计算样本均值
x_bias = np.sum(x_i-x_mean)  # 计算偏差和
print("样本均值为：",x_mean)
print("偏差和为：",x_bias)

## （2）从总体分布中抽取样本，假设总体分布为N（0，1）
x_i = np.random.randn(10)  # 从总体分布N(0,1)中抽取10个样本
x_mean = np.mean(x_i) # 计算样本均值
x_bias = np.sum(x_i-x_mean)  # 计算偏差和
print("样本均值为：",x_mean)
print("偏差和为：",x_bias)

def GetSampleDist(n, X):
    x_mean_list = []
    for i in range(n):
        x_i = np.random.choice(X, 5, replace=False)
        x_mean = np.mean(x_i)
        x_mean_list.append(x_mean)
    plt.hist(x_mean_list,color='blue',alpha=0.6,density=True)
    plt.xlabel("x")
    plt.ylabel("频率 p")
    plt.title("n="+str(n))
    plt.show()

X = np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20])   # 假设总体为X, size=20
GetSampleDist(5, X)
GetSampleDist(10, X)
GetSampleDist(20, X)
GetSampleDist(100, X)
GetSampleDist(1000, X)
GetSampleDist(10000, X)
GetSampleDist(100000, X)

从总体/总体的分布中抽取样本并计算样本方差与样本标准差

## （1）从总体中抽取样本
X = np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20])   # 假设总体为X
x_i = np.random.choice(X, 10, replace=False)  # 从总体X中抽取10个样本
x_sn2 = np.var(x_i,ddof=0) #样本方差
x_s2 = np.var(x_i,ddof=1) # 无偏样本方差
x_sn = np.std(x_i,ddof=0) # 样本标准差
x_s = np.std(x_i,ddof=1) # 无偏样本标准差
print("样本方差sn^2为：",x_sn2)
print("样本方差s^2为：",x_s2)
print("样本标准差sn^2为：",x_sn)
print("样本标准差s^2为：",x_s)

## （2）从总体分布中抽取样本，假设总体分布为N（0，1）
x_i = np.random.randn(10)  # 从总体分布N(0,1)中抽取10个样本
x_sn2 = np.var(x_i,ddof=0) #样本方差
x_s2 = np.var(x_i,ddof=1) # 无偏样本方差
x_sn = np.std(x_i,ddof=0) # 样本标准差
x_s = np.std(x_i,ddof=1) # 无偏样本标准差
print("样本方差sn^2为：",x_sn2)
print("样本方差s^2为：",x_s2)
print("样本标准差sn^2为：",x_sn)
print("样本标准差s^2为：",x_s)

x_mean_list=[]
for i in range(10000):
    x_i=np.random.randn(10)
    x_mean_list.append(np.mean(x_i))
print("标准正态分布的均值和方差为：",0,1)
print("1000个样本均值的样本均值为：",np.mean(x_mean_list))
print("1000个样本均值的样本方为：",np.var(x_mean_list,ddof=1))

最小次序统计量及其分布

from scipy.stats import rv_discrete  # 自定义离散分布
x_k = np.arange(3)+1
p_k = np.array([1/3]*3)
X = rv_discrete(name='min', values=(x_k, p_k))
def Get_Min_Dist(n, X):
    Min_list = []
    for i in range(n):
        min_value = np.min(X.rvs(size=3))
        Min_list.append(min_value)
    # 统计每个值出现的次数
    xk_count = []
    for v in x_k:
        xk_count.append(np.sum(Min_list==v))
    # 画图
    plt.bar(x_k,np.array(xk_count)/n,color='blue',alpha=0.6)
    plt.xlabel("x")
    plt.ylabel("频率 p")
    plt.xlim(0,4)
    plt.title("n="+str(n))
    plt.show()

Get_Min_Dist(5, X)
Get_Min_Dist(10, X)
Get_Min_Dist(20, X)
Get_Min_Dist(100, X)
Get_Min_Dist(1000, X)
Get_Min_Dist(10000, X)

最大次序统计量及其分布

from scipy.stats import rv_discrete  # 自定义离散分布
x_k = np.arange(3)+1
p_k = np.array([1/3]*3)
X = rv_discrete(name='min', values=(x_k, p_k))
def Get_Max_Dist(n, X):
    Max_list = []
    for i in range(n):
        max_value = np.max(X.rvs(size=3))
        Max_list.append(max_value)
    # 统计每个值出现的次数
    xk_count = []
    for v in x_k:
        xk_count.append(np.sum(Max_list==v))
    # 画图
    plt.bar(x_k,np.array(xk_count)/n,color='blue',alpha=0.6)
    plt.xlabel("x")
    plt.ylabel("频率 p")
    plt.xlim(0,4)
    plt.title("n="+str(n))
    plt.show()

Get_Max_Dist(5, X)
Get_Max_Dist(10, X)
Get_Max_Dist(20, X)
Get_Max_Dist(100, X)
Get_Max_Dist(1000, X)
Get_Max_Dist(10000, X)

从总体/总体的分布中抽取样本并计算样本中位数

## （1）从总体中抽取样本
X = np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20])   # 假设总体为X
x_i = np.random.choice(X, 10, replace=False)  # 从总体X中抽取10个样本
x_mid = np.median(x_i) # 计算样本中位数
print("样本中位数为：",x_mid)

## （2）从总体分布中抽取样本，假设总体分布为N（0，1）
x_i = np.random.randn(10)  # 从总体分布N(0,1)中抽取10个样本
x_mid = np.median(x_i) # 计算样本中位数
print("样本中位数为：",x_mid)

从总体/总体的分布中抽取样本并计算样本分位数

## （1）从总体中抽取样本
X = np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20])   # 假设总体为X
x_i = np.random.choice(X, 10, replace=False)  # 从总体X中抽取10个样本
x_low = np.percentile(x_i,25) # 计算样本下四分位数
x_high = np.percentile(x_i,75) # 计算样本上四分位数
print("样本下四分位数为：",x_low)
print("样本上四分位数为：",x_high)

def GetMidDist(n, X):
    x_mid_list = []
    for i in range(n):
        x_i = np.random.choice(X, 10, replace=False)
        x_mid = np.median(x_i)
        x_mid_list.append(x_mid)
    plt.hist(x_mid_list,color='blue',alpha=0.6,density=True)
    plt.xlabel("x")
    plt.ylabel("频率 p")
    plt.title("n="+str(n))
    plt.show()

X = np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20])   # 假设总体为X, size=20
GetMidDist(5, X)
GetMidDist(10, X)
GetMidDist(20, X)
GetMidDist(100, X)
GetMidDist(1000, X)
GetMidDist(10000, X)
GetMidDist(100000, X)

用python画出不同自由度下的卡方分布的密度函数图chi2.pdf(x,df=4)、chi2.pdf(x,df=6)、chi2.pdf(x,df=10)

# 使用scipy计算pdf画图(非自定义函数)
from scipy.stats import chi2
x = np.linspace(0.01,30,10000)  
plt.plot(x, chi2.pdf(x,df=4),'r-', lw=5, alpha=0.6, label='chi2(4)',c='red')
plt.plot(x, chi2.pdf(x,df=6),'r-', lw=5, alpha=0.6, label='chi2(6)',c='blue')
plt.plot(x, chi2.pdf(x,df=10),'r-', lw=5, alpha=0.6, label='chi2(10)',c='orange')
plt.xlabel("X")
plt.ylabel("p (x)")
plt.legend()
plt.show()

使用卡方分布的定义演示卡方分布

from scipy.stats import norm
n = 10
chi2_list = []
for i in range(100000):
    x_i = norm.rvs(loc=0,scale=1,size=10)
    chi2_T = np.sum(np.square(x_i))
    chi2_list.append(chi2_T)
sns.distplot(chi2_list,color='blue')
plt.xlabel("x")
plt.ylabel("频率 p")
plt.title("n="+str(n))
plt.show()

卡方分布的期望值等于自由度，方差等于两倍自由度

# 假设总体为N（0，1），抽样的样本容量为n=11,抽样的次数为N
from scipy.stats import norm
def S2_Chi2(N):
    mu,sig = 0,1
    n = 11
    T_list = []
    for i in range(N):
        x_i = norm.rvs(loc=mu,scale=sig,size=n)  # 正态分布总体抽样
        T = (n-1)*np.var(x_i,ddof=1)/sig**2  # 构造卡方统计量
        T_list.append(T)
    sns.distplot(T_list,color='blue')
    plt.xlabel("x")
    plt.ylabel("频率 p")
    plt.title("chi2(10)")
    plt.show()

S2_Chi2(100000)

F统计量与F分布

F是自由度为m与n的F分布

使用scipy与matplotlib绘制不同的m，n下的F分布的密度函数

from scipy.stats import f
x = np.linspace(0.01,5,10000)  
plt.plot(x, f.pdf(x,4,4000),'r-', lw=5, alpha=0.6, label='F(4,4000)',c='red')
plt.plot(x, f.pdf(x,4,10),'r-', lw=5, alpha=0.6, label='F(4,10)',c='blue')
plt.plot(x, f.pdf(x,4,4),'r-', lw=5, alpha=0.6, label='F(4,4)',c='orange')
plt.plot(x, f.pdf(x,4,1),'r-', lw=5, alpha=0.6, label='F(4,1)',c='yellow')
plt.xlabel("X")
plt.ylabel("p (x)")
plt.legend()
plt.show()

使用F统计量的定义演示

from scipy.stats import norm
m,n = 4,4000
F_list = []
for i in range(100000):
    chi2_m_sample = np.sum(np.square(norm.rvs(loc=0,scale=1,size=m))) # 卡方m统计量
    chi2_n_sample = np.sum(np.square(norm.rvs(loc=0,scale=1,size=n))) # 卡方n统计量
    F_T = (chi2_m_sample/m) / (chi2_n_sample/n)  # # F（m，n）统计量
    F_list.append(F_T)
sns.distplot(F_list,color='blue')
plt.xlabel("x")
plt.ylabel("频率 p")
plt.title("F(4,4000)")
plt.show()

F分布的密度函数是一个只取正值的偏态分布

# x1,x2,...,xn ~ N(0,1), y1,y2,...,yn ~ N(0,4)
from scipy.stats import norm
F_list = []
for i in range(100):
     x_mu, x_sigma2 = 0, 1
     y_mu, y_sigma2 = 0, 4
     norm_xi = norm.rvs(loc=x_mu,scale=x_sigma2,size=5)
     norm_yi = norm.rvs(loc=y_mu,scale=y_sigma2,size=5)
     sx_2 = np.var(norm_xi,ddof=1)
     sy_2 = np.var(norm_yi,ddof=1)
     F_T = (sx_2/x_sigma2) / (sy_2/y_sigma2)
     F_list.append(F_list)
sns.distplot(F_list,color='red')
plt.xlabel("x")
plt.ylabel("频率 p")
plt.title("F(4,4)")
plt.show()

t分布及其统计量

使用scipy与matplotlib绘制不同的n下的t分布的密度函数

from scipy.stats import t
from scipy.stats import norm 
x = np.linspace(-6,6,10000)  
plt.plot(x, t.pdf(x,4),'--', lw=5, alpha=0.6, label='t (4)',c='red')
plt.plot(x, norm.pdf(x,loc=0,scale=1),'r-', lw=5, alpha=0.6, label='N (0,1)',c='yellow')
plt.plot(x, t.pdf(x,100),'--', lw=5, alpha=0.6, label='t (100)',c='blue')
plt.xlabel("X")
plt.ylabel("p (x)")
plt.legend()
plt.show()

当自由度较大的时候，t分布近似于标准正态分布

from scipy.stats import norm 
from scipy.stats import t 
t_list = []
for i in range(300000):
    mu,sigma2 = 0,1
    x_i = norm.rvs(loc=mu, scale=sigma2, size=5)
    x_mean = np.mean(x_i)
    x_s = np.std(x_i,ddof=1)
    t_T = np.sqrt(4)*(x_mean-mu) / x_s
    t_list.append(t_T)
sns.distplot(t_list,color='blue',label='t (4)')
x = np.linspace(-6,6,10000)  
plt.plot(x, norm.pdf(x,loc=0,scale=1),'r-', lw=5, alpha=0.6, label='N (0,1)',c='yellow')
plt.xlabel("x")
plt.ylabel("频率 p")
plt.title("t (4)")
plt.xlim(-6,6)
plt.legend()
plt.show()

参数估计——点估计

矩估计

原点矩和中心距

X的k阶原点矩

X的k阶中心矩

数学期望是随机变量的1阶原点矩，方差是随机变量的二阶中心矩，原点矩刻画了随机变量X偏离原点（0，0）的程度，中心矩刻画了随机变量偏离中心的程度
假设总体是标准正态分布，求3阶原点矩和中心矩

from scipy.stats import norm
x_i = norm.rvs(loc=0, scale=1, size=10000)
a3 = np.mean(np.power(x_i,3))
b3 = np.mean(np.power((x_i-np.mean(x_i)), 3))
print("3阶原点矩：",a3)
print("3阶中心矩：",b3)

矩估计
用样本矩替代总体矩

# 假设真实值lambda = 5
from scipy.stats import expon
real_lmd = 5
x_i = np.random.exponential(scale=1/real_lmd, size=1000)
print("矩估计为：",1/np.mean(x_i))

极大似然估计

# 使用sympy演示极大似然估计的案例
from sympy import *
p = Symbol('p')  #定义总体参数
P_p = p**7*(1-p)**3  # 定义似然函数
lnP_p = ln(P_p) # 化简为对数似然
d_ln_P = diff(lnP_p, p) # 求导函数
p_hat = solve(d_ln_P, p) # 导函数为0
print("p的极大似然估计为：",p_hat)

指数函数中lamda的极大似然估计

from sympy.abc import lamda
x_1,x_2,x_3,x_4,x_5 = symbols('x_1:6') # 定义多个样本变量
x_1,x_2,x_3,x_4,x_5 = 2, 3, 0.5, 5, 2
f_lmd = lamda*E**(-lamda*x_1) * lamda*E**(-lamda*x_2) * lamda*E**(-lamda*x_3) * lamda*E**(-lamda*x_4) * lamda*E**(-lamda*x_5) # 定义似然函数
ln_f_lmd = ln(f_lmd) # 定义对数似然函数
d_ln_f = diff(ln_f_lmd, lamda) # 求导
lmd_hat = solve(d_ln_f, lamda) # 导数为0
print("指数分布参数lamda的极大似然估计值为：",lmd_hat)
print("指数分布参数lamda的极大似然公式求解为：", 5 / (x_1+x_2+x_3+x_4+x_5))

无偏性与有效性（方差最小）

区间估计

如何通过样本构造一个置信区间？？？枢轴量在很多情况下不怎么快速可行，因此可以使用bootstrap方法去构造置信区间

bootstrap就是在样本中重抽样，每抽样一次计算一次统计量T
，这样就可以构造一个抽样分布，取抽样分布的a/2分位数点和1-a/2分位数点作为区间左右边界就可

# 使用bootstrap方法计算N(0,1)的mu的置信区间：样本量为1000，重抽样样本量为500, 重抽样的次数为100000次
T_list = []
N, N_re = 1000, 500
total_times = 10000
alpha = 0.05
x_i = np.random.randn(N)  # 抽样1000个
for i in range(total_times):
    x_re = np.random.choice(x_i, N_re, replace=True) # 从样本中重抽样
    T = np.mean(x_re)
    T_list.append(T)
left = np.percentile(np.array(T_list), 100*alpha/2)
right = np.percentile(np.array(T_list), 100*(1-alpha/2))
print("正态总体的mu的置信区间为：["+str(left)+", "+str(right)+"]")

置信区间的置信水平是什么意思？？？

# 探索置信度1-alpha的含义：
def get_confident_interval(x_i):
    T_list = []
    N, N_re = 1000, 500
    total_times = 10000
    alpha = 0.05
    x_i = np.random.randn(N) 
    for i in range(total_times):
        x_re = np.random.choice(x_i, N_re, replace=True) # 从样本中重抽样
        T = np.mean(x_re)
        T_list.append(T)
    left = np.percentile(np.array(T_list), 100*alpha/2)
    right = np.percentile(np.array(T_list), 100*(1-alpha/2))
    return {'left':left, 'right':right}

left_right_list = []  #100个置信区间的列表
for i in range(100):
    x_i = np.random.randn(1000) #每次抽样1000个样本
    T_i = get_confident_interval(x_i)
    left_right_list.append(T_i)
for i in range(len(left_right_list)):
    plt.vlines(x=i+1, ymin=left_right_list[i]['left'], ymax=left_right_list[i]['right'])
    plt.scatter(np.array([i+1]*2),np.array([left_right_list[i]['left'],left_right_list[i]['right']]),color='blue')
plt.axhline(y=0, ls='--', color='b')
plt.show()

假设检验之基本思想

步骤（建立原假设与备择假设——选择统计量并给出拒绝域的形式——选择显著性水平a——给出拒绝域W）

# python模拟以上的假设检验：
from scipy.stats import norm
x_i = np.array([1.0, 0.85, 0.90, 0.98, 0.96, 1.0, 0.97, 0.98, 0.98, 0.98]) # 从当天生产的肉类总体中抽样10个
x_mean = np.mean(x_i) # 构造样本均值统计量
c = norm(loc=1, scale=np.sqrt(1/10)).ppf(0.05)  # 总体N(theta, 1)，则样本均值的抽样分布为N(theta, 1/n)，n为样本量
print("样本均值为：",x_mean)
print("拒绝域为：（-inf， "+str(c)+"]")
print("x_mean位于接受域，接受原假设！")

假设检验之正态总体参数的假设检验

方差已知时的u检验

# 使用python模拟u检验：
from scipy.stats import norm
mu_0, sigma= 65, 0.4
x_i = np.array([60, 72, 67, 89, 90, 100, 67, 76, 87, 78])  #样本
u_T = (np.mean(x_i) - mu) / (sigma / np.sqrt(len(x_i))) # 计算样本统计量
left = norm(loc=0, scale=1).ppf(0.025)
right = norm(loc=0, scale=1).ppf(0.975)
print("拒绝域为：["+str(left)+","+str(right)+"]")
print("u统计量为：",u_T)
print("u_T位于拒绝域，拒绝原假设！")

方差未知时的t检验

# 使用python模拟单样本t检验
from scipy.stats import t
mu_0 = 250
x_i = np.array([248.8, 249.2, 250.7, 251.2, 248.0, 253.0, 248.9, 250.2, 251.2, 249.2]) # 样本
t_T = (np.mean(x_i)-mu_0)/(np.std(x_i, ddof=1)/np.sqrt(len(x_i)-1)) # 样本统计量
right = t(len(x_i)-1).ppf(0.05)
print("拒绝域为：(-inf, "+str(right)+"]")
print("t统计量为：",t_T)
print("t_T位于接受域，不能拒绝原假设！")

# 使用scipy进行单样本t检验
from scipy.stats import ttest_1samp
x_i = np.array([248.8, 249.2, 250.7, 251.2, 248.0, 253.0, 248.9, 250.2, 251.2, 249.2]) # 样本
t_T, p_value = ttest_1samp(x_i, popmean=250, alternative="less") # popmean:总体的期望， alternative={‘two-sided’（双侧）, ‘less’（左侧）, ‘greater’（右侧）} ； 返回：t统计量和p值
print("p = ", p_value)
print("由于p>0.05，不能拒绝原假设！")

两个正态总体的均值差检验(方差已知)

# 使用python进行两个正态总体的均值差检验
from scipy.stats import norm 
sigma2_1 = sigma2_2 = 0.4
x_i = np.array([2.2, 3.1, 2.5, 2.8, 3.0, 5.0, 4.3, 1.8, 0.2, 1.0]) # 修复前的样本
y_i = np.array([3.2, 2.8, 2.2, 1.4, 4.4, 3.4, 6.4, 4.0, 3.9, 2.1]) # 修复后的样本
m, n = len(x_i), len(y_i)  # 样本数
u_T = (np.mean(x_i) - np.mean(y_i)) / np.sqrt(sigma2_1/m + sigma2_2/n)  #计算检验统计量
left = norm(loc=0, scale=1).ppf(0.025)  #拒绝域左边界
right = norm(loc=0, scale=1).ppf(0.975) # 拒绝域右边界
print("样本统计量为：", u_T)
print("拒绝域为：(-inf, "+str(left)+"]与["+str(right)+", +inf)")
print("样本统计量位于拒绝域，拒绝原假设！")

方差未知时的两样本t检验

# 使用python的scipy进行方差未知时，两个正态总体的均值差检验
from scipy.stats import ttest_ind
x_i = np.array([2.2, 3.1, 2.5, 2.8, 3.0, 5.0, 4.3, 1.8, 0.2, 1.0]) # 修复前的样本
y_i = np.array([3.2, 2.8, 2.2, 1.4, 4.4, 3.4, 6.4, 4.0, 3.9, 2.1]) # 修复后的样本
m, n = len(x_i), len(y_i)  # 样本数
t_T, p_value = ttest_ind(x_i, y_i, equal_var=False, alternative='two-sided') # equal_var=True默认同方差， alternative={‘two-sided’（双侧）, ‘less’（左侧）, ‘greater’（右侧）} ； 返回：t统计量和p值
print("p值为：",p_value)
print("p>0.05, 则不能拒绝原假设！")

单个正态总体方差的检验

# 使用python进行单个正态总体的方差检验
from scipy.stats import chi2
sigma0_2 = 2.0
x_i = np.array([23.3, 24.7, 23.1, 19.4, 18.9, 20.3, 24.9, 30.2, 19.0, 20.9]) # 样本
n = len(x_i) # 样本量
chi_T = (n-1)*np.std(x_i, ddof=1) / sigma0_2 #检验统计量
W_right = chi2(df = n-1).ppf(0.95) #拒绝域的边界
print("检验统计量为：",chi_T)
print("拒绝域为：["+str(W_right)+", +inf)")
print("统计量位于接受域，不能拒绝原假设！")

两个正态总体方差比的 F检验

# 使用python做两个正态总体方差比的 $F$ 检验：
from scipy.stats import f
x_i = np.array([23.3, 24.7, 23.1, 19.4, 18.9, 20.3, 24.9, 30.2, 19.0, 20.9]) # 样本X
y_i = np.array([34.1, 43.2, 31.1, 12.4, 36.9, 21.3, 39.9, 20.2, 39.0, 12.2]) # 样本Y
m,n = len(x_i), len(y_i) # 样本量
F_T = np.std(x_i, ddof=1) / np.std(y_i, ddof=1) #检验统计量
W_left = f(dfn=m-1, dfd=n-1).ppf(0.05)
print("统计量为：", F_T)
print("拒绝域为：(0, "+str(W_left)+"]")
print("统计量位于拒绝域，拒绝原假设！")

假设检验之似然比检验与Bootstrap方法

似然比检验

似然比检验的统计量并没有一个统一的精确分布，但是2倍对数似然比渐近服从卡方分布，卡方分布的自由度为独立参数的个数

# 使用Bootstrap做假设检验
mu_0 = 250
x_i = np.array([248.8, 249.2, 250.7, 251.2, 248.0, 253.0, 248.9, 250.2, 251.2, 249.2]) # 样本
x_bar_T = np.mean(x_i)  #计算检验统计量
## Bootstrap重抽样构造统计量的抽样分布过程
T_list = []
N = 10000  # 重抽样次数
for i in range(N):
    x_i_resample = np.random.choice(x_i, len(x_i), replace=True)
    x_bar_resample = np.mean(x_i_resample)
    T_list.append(x_bar_resample)
right = np.percentile(T_list, 5)
print("检验统计量为：",x_bar_T)
print("拒绝域为：(-inf, "+str(right)+"]")
print("检验统计量位于接受域，不能拒绝原假设！")

# t检验解决以上问题
from scipy.stats import t
mu_0 = 250
x_i = np.array([248.8, 249.2, 250.7, 251.2, 248.0, 253.0, 248.9, 250.2, 251.2, 249.2]) # 样本
t_T = (np.mean(x_i)-mu_0)/(np.std(x_i, ddof=1)/np.sqrt(len(x_i)-1)) # 样本统计量
right = t(len(x_i)-1).ppf(0.05)
print("拒绝域为：(-inf, "+str(right)+"]")
print("t统计量为：",t_T)
print("t_T位于接受域，不能拒绝原假设！")