概率论与数理统计
一、描述性统计和统计图
1.用Pandas来计算统计量
使用 pandas的describe方法计算相关统计量，并计算身高和体重的偏度，峰度，样本的25%,50%,90%分位数
数据如上图所示
from numpy import reshape,c_
import pandas as pd
df = pd.read_excel('F:/hellopython/Pdata4_6_1.xlsx',header=None)
a=df.values
h=a[:,::2];#身高
w=a[:,1::2]#体重
h=reshape(h,(-1,1)) #弄成一列
w=reshape(w,(-1,1))
df2=pd.DataFrame(c_[h,w],columns=['身高','体重'])
print(df2.describe())
print('偏度为',df2.skew())
print('峰度为',df2.kurt())
print('分位数为',df2.quantile(0.9))
2.统计图
2.1频数图及直方图
计算频数并且画直方图的命令为：
hist(x,bins=None,range=None,density=None,weights=None,cumulative=False,bottom=None,histtype=‘bar’ ,align='mid ’ ,orientation='vertical ’ , rwidth=None，logFalse,color=None,label=None,stacked=False)
他将区间[min(x),max(x)]等分为bins份，统计在每个左闭右开小区间（最后一个区间为闭区间）上的数据出现的频数并画直方图。其中一些参数含义如下：
接上题画出身高和体重的直方图，并统计从最小体重到最大体重，等间距分为6个小区间，数据出现在每个小区间的频数
import numpy as np
import matplotlib.pyplot as plt
df = pd.read_excel('F:/hellopython/Pdata4_6_1.xlsx',header=None)
h=a[:,::2]; w=a[:,1::2]
h=np.reshape(h,(-1,1)); w=np.reshape(w,(-1,1))
plt.rc('font',size=16); plt.rc('font',family="SimHei")
plt.subplot(121); plt.xlabel("身高"); plt.hist(h,10) #只画直方图不返回频数表
plt.subplot(122); ps=plt.hist(w,6) #画图并返回频数表ps
plt.xlabel("体重"); print("体重的频数表为：", ps)
plt.savefig("figure4_8.png", dpi=500); plt.show()
体重的频数表为： (array([ 9., 13., 27., 31., 11., 9.]), array([47., 52., 57., 62., 67., 72., 77.]), <BarContainer object of 6 artists>)
2.2箱线图
下面分别给出了25个男子和25个女子的肺活量数据（以L计，数据已经排过序）（略），请画出线相图
import numpy as np
import matplotlib.pyplot as plt
a=np.loadtxt("Pdata4_9.txt") #读入两行的数据
b=a.T #转置成两列的数据
plt.rc('font',size=16); plt.rc('font',family='SimHei')
plt.boxplot(b,labels=['女子','男子'])
plt.savefig('figure4_9.png', dpi=500); plt.show()
箱线图特别适用于比较两个或两个以上数据集的性质,为此,将几个数据集的箱线图画在同一个图形界面上.例如,在图4.4中可以明显地看到男子的肺活量要比女子的肺活量大,男子的肺活量较女子的肺活量分散.
箱线图异常值检测
在数据集中某一个观察值不寻常地大于或小于该数集中的其他数据,称为疑似异常值.疑似异常值的存在,会对随后的计算结果产生不适当的影响.检查疑似异常值并加以适当的处理是十分必要的.
第一四分位数Q1与第三四分位数Q3之间的距离:Q3-Q1记为_IOR,称为四分位数间距．若数据小于Q1-1.5IQR或大于Q3＋ 1.5IQR,就认为它是疑似异常值.
看看身高和体重的箱线图
import numpy as np
import matplotlib.pyplot as plt
a=np.loadtxt("Pdata4_6_2.txt")
h=a[:,::2]; w=a[:,1::2]
h=np.reshape(h,(-1,1)); w=np.reshape(w,(-1,1))
hw=np.hstack((h,w)); plt.rc('font',size=16)
plt.rc('font',family='SimHei')
plt.boxplot(hw,labels=['身高','体重'])
plt.savefig("figure4_10.png",dpi=500); plt.show()
2.3经验分布函数
画出体重的经验分布函数图形
import numpy as np
import matplotlib.pyplot as plt
a=np.loadtxt("Pdata4_6_2.txt")
w=a[:,1::2]; w=np.reshape(w,(-1,1)); plt.rc('font',size=16)
h=plt.hist(w,20,density=True, histtype='step', cumulative=True)
print(h); plt.grid()
plt.savefig("figure4_11.png",dpi=500); plt.show()
2.4 QQ图
QQ图可以观测两组数据是否属于同一分布，或者观测一组数据是否属于正态分布
对于例4.6中的身高数据，如果他们来自于正态分布，求该矩估计正态分布的参数，试画出他们的Q-Q图，判断拟合效果。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, probplot
a=np.loadtxt("Pdata4_6_2.txt")
h=a[:,::2]; h=h.flatten()
mu=np.mean(h); s=np.std(h); print([mu,s])#求样本均值和方差
sh=np.sort(h) #按从小到大排序
n=len(sh); xi=(np.arange(1,n+1)-1/2)/n
yi=norm.ppf(xi,mu,s)
plt.rc('font',size=16);plt.rc('font',family='SimHei')
plt.rc('axes',unicode_minus=False) #用来正常显示负号
plt.subplot(121); plt.plot(yi, sh, 'o', label='QQ图');
plt.plot([155,185],[155,185],'r-',label='参照直线')
plt.legend(); plt.subplot(122)
res = probplot(h,plot=plt)
plt.savefig("figure4_12.png",dpi=500); plt.show()
二、参数估计和检验假设
1.参数估计
1.1极大似然估计
假定例4.6中学生的身高服从正态分布，求总体均值和标准差的极大似然估计。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
a=np.loadtxt("Pdata4_6_2.txt")
h=a[:,::2]; h=h.flatten()
mu=np.mean(h); s=np.std(h);
print("样本均值和标准差为：",[mu,s])
print("极大似然估计值为：", norm.fit(h))
1.2区间估计
直接讲例子把，有一大批糖果，现从中选择16袋，称得重量（以g计）如下（略），设袋装糖果的重量近似的服从正态分布，试求总体均值 mu 的置信水平为0.95的置信区间。
from numpy import array, sqrt
from scipy.stats import t
a=array([506, 508, 499, 503, 504, 510, 497, 512,
514, 505, 493, 496, 506, 502, 509, 496])
# ddof取值为1时，标准偏差除的是(N-1)；NumPy中的std计算默认是除以N
mu=a.mean(); s=a.std(ddof=1) #计算均值和标准差
print(mu, s); alpha=0.05; n=len(a)
val=(mu-s/sqrt(n)*t.ppf(1-alpha/2,n-1),mu+s/sqrt(n)*t.ppf(1-alpha/2,n-1))
print("置信区间为：",val)
或者直接调用库
import numpy as np
import scipy.stats as ss
from scipy import stats
a=np.array([506, 508, 499, 503, 504, 510, 497, 512,
514, 505, 493, 496, 506, 502, 509, 496])
alpha=0.95; df=len(a)-1
ci=ss.t.interval(alpha,df,loc=a.mean(),scale=ss.sem(a))
print("置信区间为：",ci)
2.参数假设检验（知道总体服从什么类型的分布）
2.1. 单个总体均值的假设检验
（1）正态总体标准差 sigma 已知的Z检验法
例：某车间用一台包装机包装糖果.包得的袋装糖重是一个随机变量,它服从正态分布.当机器正常时,其均值为0.5kg,标准差为0.015kg.某日开工后为检验包装机是否正常,随机地抽取它所包装的糖9袋,称得净重(kg)为
0.497 0.506 0.518 0.524 0.498 0.511 0.520 0.515 0.512问机器是否正常?
import numpy as np
from statsmodels.stats.weightstats import ztest
sigma=0.015
a=np.array([0.497, 0.506, 0.518, 0.524, 0.498, 0.511, 0.520, 0.515, 0.512])
tstat1, pvalue=ztest(a,value=0.5) #计算T统计量的观测值及p值
tstat2=tstat1*a.std(ddof=1)/sigma #转换为Z统计量的观测值
print('t值为：',round(tstat1,4))
print('z值为：',round(tstat2,4)); print('p值为:',round(pvalue,4))
（2）正态总体标准差 \sigma 未知的 t 检验法
例： 某批矿砂的5个样品中的镍含量（%),经测定为3.25 3.27 3.24 3.26 3.24
设测定值总体服从正态分布,但参数均未知,问在α=0.01下能否接受假设:这批矿砂的镍含量的均值为3.25.
import numpy as np
from statsmodels.stats.weightstats import ztest
a=np.array([3.25, 3.27, 3.24, 3.26, 3.24])
tstat, pvalue=ztest(a,value=3.25)
print('检验统计量为：',tstat); print('p值为:',pvalue)
2.对两个样本的检验（方差相等
两位作家作品中单词统计数据
马克·吐温 0.225 0.262 0.217 0.240 0.230 0.229 0.235 0.217
斯诺特格拉斯 0.209 0.205 0.196 0.210 0.202 0.207 0.224 0.223 0.220 0.201
设两组数据分别来自正态总体,且两总体方差相等,但参数均未知.两样本相互独立．问两位作家所写的小品文中包含由3个字母组成的单词的比例是否有显著的差异(取α =0.05)?
import numpy as np
from statsmodels.stats.weightstats import ttest_ind
a=np.array([0.225, 0.262, 0.217, 0.240, 0.230, 0.229, 0.235, 0.217])
b=np.array([0.209, 0.205, 0.196, 0.210, 0.202, 0.207,
0.224, 0.223, 0.220, 0.201])
tstat, pvalue, df=ttest_ind(a, b, value=0)
print('检验统计量为：',tstat); print('p值为:',pvalue)
print('自由度为：',df)
3.非参数假设检验（不知道总体服从什么类型的分布）
1. 分布拟合检验（ ka^2 检验法）
根据某市公路交通部门某年上半年交通事故记录,统计得星期一至星期日发生交通事故的次数如表4.8所示.试检验交通事故的发生次数是否服从离散均匀分布,即交通事故的发生与星期几无关.
import numpy as np
import scipy.stats as ss
bins=np.arange(1,8)
mi=np.array([36, 23, 29, 31, 34, 60, 25])
n=mi.sum(); p=np.ones(7)/7
cha=(mi-n*p)**2/(n*p); st=cha.sum()
bd=ss.chi2.ppf(0.95,len(bins)-1) #计算上alpha分位数
print("统计量为：{}，临界值为：{}".format(st,bd))
是否符合正态分布
某车间生产滚珠,随机地抽出50粒,测符它的且江(平: )为
15.0 15.8 15.2 15.1 15.9 14.7 14.8 15.5 15.6 15.3
15.1 15.3 15.0 15.6 15.7 14.8 14.5 14.2 14.9 14.9
15.2 15.0 15.3 15.6 15.1 14.9 14.2 14.6 15.8 15.2
15.9 15.2 15.0 14.9 14.8 14.5 15.1 15.5 15.5 15.1
15.1 15.0 15.3 14.7 14.5 15.5 15.0 14.7 14.6 14.2
经过计算知样本均值 = 15.0780,样本标准差s=0.4282,试问滚珠直径是否服从正态分布
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
n=50; k=8 #初始小区间划分的个数
a=np.loadtxt("Pdata4_20.txt")
a=a.flatten(); mu=a.mean(); s=a.std()
print("均值为：", mu); print("标准差为：", s)
print("最大值为：",a.max()); print("最小值为：",a.min())
bins=np.array([14.2, 14.625, 14.8375, 15.05, 15.2625, 15.475, 15.9])
h=plt.hist(a,bins)
f=h[0]; x=h[1] #提取各个小区间的频数和小区间端点的取值
print("各区间的频数为：",f,"\n小区间端点值为：",x)
p=ss.norm.cdf(x, mu, s) #计算各个分点分布函数的取值
dp=np.diff(p) #计算各小区间取值的理论概率
dp[0]=ss.norm.cdf(x[1],mu,s) #修改第一个区间的概率值
dp[-1]=1-ss.norm.cdf(x[-2],mu,s) #修改最后一个区间的概率值
print("各小区取值的理论概率为：",dp)
st=sum(f**2/(n*dp))-n #计算卡方统计量的值
bd=ss.chi2.ppf(0.95,k-5) #计算上alpha分位数
print("统计量为：{}，临界值为：{}".format(st,bd))
2. Kolmogorov-Smirnov检验
检验学生的体重是否服从正态分布
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
a=np.loadtxt("Pdata4_6_2.txt")
w=a[:,1::2]; w=w.flatten()
mu=w.mean(); s=w.std(ddof=1) #计算样本均值和标准差
print("均值和标准差分别为：", (mu, s))
statVal, pVal=ss.kstest(w,'norm',(mu,s))
print("统计量和P值分别为：", [statVal, pVal])
4.方差分析
方差分析是用于两个及两个以上总体均值差别的显著性检验。例如，生产某种产品，为了使生产过程稳定，达到优质、高产，需要对影响产品质量的因素进行分析，找出有显著影响的那些因素，方差分析就是鉴别各因素效应的一种有效的统计方法
1.单因素方法分析
若在一项试验中，只考虑一个因素A的变化，其他因素保持不变，称这种试验为单因素试验。具体做法是：A取几个水平，在每个水平上作若干个试验，试验过程中除A外其他影响指标的因素都保持不变（只有随机因素存在），我们的任务是从试验结果推断，因素A对指标有无显著影响，即当A取不同水平时指标有无显著差别。A取某个水平下的指标可视为随机变量，判断A取不同水平时指标有无显著差异，相当于检验若干总体的均值是否相等
用4种工艺A1, A2, A3,A4 生产灯泡,从各种工艺制成的灯泡中各抽出了若干个测量其寿命,结果如表4.10所示，试推断这几种工艺制成的灯泡寿命是否有显著差异.
具体原理可以再去详细了解，这里直接给出Python程序
import numpy as np
import statsmodels.api as sm
y=np.array([1620, 1670, 1700, 1750, 1800, 1580, 1600, 1640, 1720,
1460, 1540, 1620, 1680, 1500, 1550, 1610])
x=np.hstack([np.ones(5), np.full(4,2), np.full(4,3), np.full(3,4)])
d= {'x':x,'y':y} #构造字典
model = sm.formula.ols("y~C(x)",d).fit() #构建模型
anovat = sm.stats.anova_lm(model) #进行单因素方差分析
print(anovat)
双因素方法分析
如果要考虑两个因素对指标的影响，就要采用双因素方差分析。它的基本思想是：对每个因素各取几个水平，然后对各因素不同水平的每个组合作一次或若干次试验，对所得数据进行方差分析。对双因素方差分析可分为无重复和等重复试验两种情况，无重复试验只需检验两因素是否分别对指标有显著影响；而对等重复试验还要进一步检验两因素是否对指标有显著的交互影响。
具体原理自行了解，这里直接给出例子和代码
表4.17给出了某种化工过程在三种浓度、四种温度水平下得率的数据。试在水平 下，检验在不同温度（因素A）、不同浓度（因素B）下的得率是否有显著差异？交互作用是否显著？
import numpy as np
import statsmodels.api as sm
y=np.array([[11, 11, 13, 10], [10, 11, 9, 12],
[9, 10, 7, 6], [7, 8, 11, 10],
[5, 13, 12, 14], [11, 14, 13, 10]]).flatten()
A=np.tile(np.arange(1,5),(6,1)).flatten()
B=np.tile(np.arange(1,4).reshape(3,1),(1,8)).flatten()
d={'x1':A,'x2':B,'y':y}
model = sm.formula.ols("y~C(x1)+C(x2)+C(x1):C(x2)",d).fit() #注意交互作用公式的写法
anovat = sm.stats.anova_lm(model) #进行双因素方差分析
print(anovat)
运行结果：
df sum_sq mean_sq F PR(>F)
C(x1) 3.0 19.125000 6.375000 1.330435 0.310404
C(x2) 2.0 40.083333 20.041667 4.182609 0.041856
C(x1):C(x2) 6.0 18.250000 3.041667 0.634783 0.701009
Residual 12.0 57.500000 4.791667 NaN NaN
其中 =0.310404 0.041856 0.701009，即认为温度影响不显著，而浓度因素有显著差异，交互作用不显著。
1677
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
