Chp6-2
In [2]: import numpy as np import pandas as pd from scipy import stats import matplotlib.pyplot as plt from scipy.stats import t np.random.seed(1234) my_data1=stats.poisson.rvs(loc=10,mu=60,size=3000) # 生成一个规定均值的泊松分布 #pd.Series(my_data1).hist().get_figure().show #print(' 第一个分布的均值是： 70,\t 统计平均是： ',my_data1.mean()) my_data2=stats.poisson.rvs(loc=10,mu=15,size=6000) #pd.Series(my_data2).hist().get_figure().show #print(' 第二个分布的均值是： 25，统计平均是： ',my_data2.mean()) my_data=np.concatenate((my_data1,my_data2)) # 一个典型的 bi-model 分布，我们以这 9000 个数作为总体 #print(' 总体的均值为: ',my_data.mean()) # 假装我们不知道 # 现在根据样本，验证 H0: 总体的均值是 47.5 # 这是单样本双边 t 检验， one-sample t-test print('空假设 H0 是：总体的均值是 47.5 \n') sample_data=np.random.choice(a=my_data,size=100) # 从中随机抽取 100 个做样本 t_statistic,p_value=stats.ttest_1samp(a=sample_data,popmean=47.5) print('从样本构造的 t 统计量 = ',t_statistic) #t 的绝对值越大， p 值越小，可参见下面的 t-分布图 print('单样本双边检验的 p = ',p_value) # 画一下样本容量为 100 时的 t 分布曲线 df=100-1 x = np.linspace(stats.t.ppf(0.00000001,df),stats.t.ppf(0.99999999,df),100) #ppf 函数是 CDF 的逆函数，用来求分位点 plt.plot(x, t.pdf(x, df)) #pdf 生成概率密度函数表 plt.plot((t_statistic,t_statistic),(-0.01,0.4),'-.r') str_legend=('t distribution','calculated t') plt.legend(str_legend) plt.show() 空假设 H0 是：总体的均值是 47.5 从样本构造的 t 统计量 = -3.7405281096559153 单样本双边检验的 p = 0.00030786517310847877
In [3]: x = np.linspace(stats.t.ppf(0.00000001,df),stats.t.ppf(0.99999999,df),100) #ppf 函数是 CDF 的逆函数，用来求分位点 plt.plot(x, t.pdf(x, df)) #pdf 生成概率密度函数表 plt.plot((t_statistic,t_statistic),(-0.01,0.4),'-.r') str_legend=('t distribution','calculated t') plt.legend(str_legend) x_025=stats.t.ppf(0.025,df) # 注意：同样选取 alpha=0.5，做双边 test 时要给左、右两侧各分配 0.25 plt.plot((x_025,x_025),(-0.01,0.4),'--g') # 这根线是双边检验的左阈值 x_975=stats.t.ppf(0.975,df) plt.plot((x_975,x_975),(-0.01,0.4),'--g') # 这根线是双边检验的右阈值 #print(x_975) rect1=plt.Rectangle((x_025,0),x_975*2,0.5,color='g',alpha=0.25) plt.gca().add_patch(rect1) rect2=plt.Rectangle((x_975,0),5,0.5,color='r',alpha=0.25) plt.gca().add_patch(rect2) rect3=plt.Rectangle((x_025,0),-5,0.5,color='r',alpha=0.25) plt.gca().add_patch(rect3) #plt.ylim(-0.02,0.5) plt.text(x_025,0.4,'accept region') # 图片上不支持中文 #plt.savefig('111.png') # 不支持中文文件名 plt.show()
In [36]: print('p = ',p_value,'\n') alpha=0.05 if p_value>=0.05: print('\033[1;31m 接受 \033[0m H0') else: print('\033[1;31m 拒绝 \033[0m H0, 即总体的均值不等于 47.5，此时错误拒绝 H0 的概率为 ',p_value,'小于显著性水平') p = 0.00030786517310847877 拒绝 H0, 即总体的均值不等于 47.5，此时错误拒绝 H0 的概率为 0.00030786517310847877 小于显著性水平 In [18]: import numpy as np import pandas as pd from scipy import stats from scipy.stats import t import matplotlib.pyplot as plt str_legend=('t distribution','alpha=0.05 threshod','alpha=0.025 threshold',\ 'alpha=0.975 threshold','real t') # 画一下样本容量为 100 时的 t 分布曲线 df=100-1 #a=0.05 # 自定义的显著性水平 x = np.linspace(stats.t.ppf(0.00000001,df),stats.t.ppf(0.99999999,df),100) #ppf 函数是 CDF 的逆函数，用来求分位点 plt.plot(x, t.pdf(x, df), alpha=1) # 注意此 alpha 是 plot 方法中的透明度参数，不是显著性参数 x_05=stats.t.ppf(0.05,df) #p(t<x_05)=0.05 #print(x_05) plt.plot((x_05,x_05),(-0.01,0.4),'--') # 这根线是单边检验的阈值 x_025=stats.t.ppf(0.025,df) # 注意：同样选取 alpha=0.5，做双边 test 时要给左、右两侧各分配 0.25 plt.plot((x_025,x_025),(-0.01,0.4),'--g') # 这根线是双边检验的左阈值 x_975=stats.t.ppf(0.975,df) plt.plot((x_975,x_975),(-0.01,0.4),'--g') # 这根线是双边检验的右阈值 t_real=stats.t.ppf(p_value/2,df) # 上一个 cell 中我们做的是双边 t-test，因此这里找分位点时， p 值要除以 2 print('从样本获得的 t 统计量是',t_real) plt.plot((x_real,x_real),(-0.01,0.4),'-.r') # 各阈值截断的曲线下面积对应 I 型错误率 plt.legend(str_legend) plt.show() 从样本获得的 t 统计量是 -3.740528109655912
In [58]: np.random.seed(1234) my_data1=stats.poisson.rvs(loc=10,mu=60,size=3000) # 生成一个规定均值的泊松分布 #pd.Series(my_data1).hist().get_figure().show #print(' 第一个分布的均值是： 70,\t 统计平均是： ',my_data1.mean()) my_data2=stats.poisson.rvs(loc=10,mu=15,size=6000) #pd.Series(my_data2).hist().get_figure().show #print(' 第二个分布的均值是： 25，统计平均是： ',my_data2.mean()) my_data=np.concatenate((my_data1,my_data2)) # 一个典型的 bi-model 分布，我们以这 9000 个数作为总体 my_sample={} for n in range(2): my_sample[n]=np.random.choice(a=my_data,size=100) # 从中随机抽取 100 个做样本 print('第',n,'组样本的均值为',my_sample[n].mean()) print('根据样本均值，能否下结论说两组样本来源的总体其均值不相等呢？ ') # 需进行 two sample t-test #H0: 两样本均值相等 第 0 组样本的均值为 39.3 第 1 组样本的均值为 41.14 根据样本均值，能否下结论说两组样本来源的总体均值不相等呢？ In [59]: # 能否下结论说两组样本来源的总体其均值不相等呢？ # 需进行 two sample t-test #H0: 两样本均值相等 alpha=0.01 # 设置显著性水平 t_statistic,p_value=stats.ttest_rel(a=my_sample[0],b=my_sample[1]) #sample size 相等的双样本均值比较 print('t = ',t_statistic) print('p = ',p_value) if p_value<=alpha: print('\033[1;31m 拒绝 \033[0m H0：两样本来源的总体均值相等 ') else: print('\033[1;31m 接受 \033[0m the H0: 两样本来源的总体均值相等') t = -0.5797156447793128 p = 0.5634233550606108 接受 the H0: 两样本来源的总体均值相等 In [1]: #Titanic 数据存活人的类别分析 # 探究性别与存活是否为关联事件 # 如果性别与存活与否无关，则存活人群中的性别比应与全体人的一致，全体是什么？ # 是船上所有人，还是所有人（即 0.5:0.5） import pandas as pd import numpy as np from scipy import stats from scipy.stats import chi2 from matplotlib import pyplot as plt titanic = pd.read_csv("C:\Python\Scripts\my_data\Titanic.csv") # 对全体做性别统计 mask1=titanic['Sex']=='male' mask2=titanic['Sex']=='female' p=np.array([sum(mask1)/(sum(mask1)+sum(mask2)),sum(mask2)/(sum(mask1)+sum(mask2))]) print(p) mask_survived=titanic['Survived']==0 # 如果反过来做，活着的人中，性别比满足 p 吗？请同学们验证 my_population=titanic.loc[mask_survived,'Sex'] print(type(my_population)) pop_size=my_population.count() print('总体大小是 ',pop_size) sample=my_population E=pop_size*p print( '预期的男、女个数是 ',E) mask1=sample=='male' mask2=sample=='female' my_set1=sample[mask1] my_set2=sample[mask2] O=np.array([len(my_set1),len(my_set2)]) print('实际的男女个数是 ',O) chi_squard,p_value=stats.chisquare(f_obs=O,f_exp=E) print(chi_squard,p_value) a=0.05 df=1 # 只有两个类别，因此自由度是 2-1=1 if p_value<=a: print('我们 \033[1;31m 拒绝 \033[0m 男性女性具有相同死亡率的假设.') else: print('我们 \033[1;31m 接受 \033[0m 男性女性具有相同死亡率的假设') # 请同学们尝试分析仓位等级、登船码头与存活率之间是否有关联 [0.64758698 0.35241302] <class 'pandas.core.series.Series'> 总体大小是 549 预期的男、女个数是 [355.52525253 193.47474747] 实际的男女个数是 [468 81] 100.96890721903955 9.343874013954211e-24 我们 拒绝 男性女性具有相同死亡率的假设.
