Scipy_常用统计函数

Scipy_常用统计函数

Weijia Li

关注公众号，一起学习~

无聊e编程

from scipy.stats import *
import numpy as np
import matplotlib.pyplot as plt

describe

参量(a, axis=0, ddof=1, bias=True, nan_policy=‘propagate’) a:输入数据 axis:默认0;none计算整个a ddof:自由度默认为1 bias:如果为False，则校正偏度和峰度计算以消除统计偏差 nan_policy:定义当输入包含nan时如何处理。可以使用以下选项 ‘propagate’:返回nan ‘raise’:抛出一个错误 ‘omit’:执行计算时忽略nan值

返回: nobs:观察数（沿轴的数据长度) minmax：数据数组的最小值和最大值。mean:数据沿轴的算术平均值。方差:沿轴,分母的数据的无偏方差是观察数减去一。偏度:基于分母等于观测次数的矩计算，即没有自由度校正。峰度:峰度被归一化，因此对于正态分布它为零。没有使用自由度。

集中趋势的描述

mean算术均数

# 抗体滴度为1:10,1:20，1:40，1:80，1:160，1:320，1:640，1:1280
# 分别对应人数f为：4,3,10,10,11,15,14,2
Didu=np.array([10,20,40,80,160,320,640,1280])
Repeat_f=np.array([4,3,10,10,11,15,14,2])
Data=Didu.repeat(Repeat_f)
Data.mean()

280.8695652173913

gmean几何平均值

gmean(Data)

150.64109410382258

median中位数

np.median(Data)

160.0

np.percentile(Data,[50])

array([160.])

np.percentile(Data,[50])

array([160.])

percentile百分位数

np.percentile(Data,[25,50,100])

array([ 80., 160., 1280.])

scoreatpercentile(Data,per=[25,50,100])

array([ 80., 160., 1280.])

离散趋势的描述

极差

np.ptp(Data)

1270

四分位数间距

iqr(Data)

240.0

方差，标准差

Data.std(),Data.var()

(281.48539298696903, 79234.02646502838)

变异系数

#标准差与平均值之比
variation(Data)

1.0021925756501993

其他统计描述

样本偏度

skew(Data)

1.4758210321910645

样本峰度

kurtosis(Data)

2.3411133773829924

平均值的标准误

sem(Data)

34.135117863346174

相关函数

Pearsonr相关系数

a=np.arange(7)
b=np.array([0,3,2,1,3,4,8])
plt.plot(a,b)

#return:皮尔逊相关系数和两尾p值

pearsonr(a,b)

(0.8067793113007157, 0.028319564882430495)

总体均数的估计

可信区间

#a为置信水平
#第一个参数为置信水平
t.interval(0.95,loc=166.95,df=9,scale=3.64/np.sqrt(10))

(164.34610086233263, 169.55389913766734)

#手工：
#假设样本均值为166.95，自由度为10-1，标准差为3.64,置信度为1-a
[a,μ,df,s]=[0.05,166.95,9,3.64]
print('均值的置信区间95%CI为:(',(μ-t.isf(a/2,df=df)*(s/np.sqrt(df+1))),',',(μ+t.isf(a/2,df=df)*(s/np.sqrt(df+1))),')')

均值的置信区间95%CI为:( 164.34610086233263 , 169.55389913766734 )

#若某市某年18岁男生身高服从μ=167.7，标准差δ=5.3的正态分布，总体中随机抽样g=100个，每次样本含量n=10人,置信度为1-0.05
#总体δ未知，小样本，服从t分布
g,n,μ,s,a=100,10,167.7,5.3,0.05
data=norm.rvs(loc=167.7,scale=5.3,size=n)
t.interval(0.95,loc=167.7,scale=sem(data),df=n-1)

(164.20057600017955, 171.19942399982043)

#实验1：
#若某市某年18岁男生身高服从μ=167.7，标准差δ=5.3的正态分布，总体中随机抽样g=100个，每次样本含量n=10人,置信度为1-0.05
#总体δ未知，小样本，服从t分布
g,n,μ,s,a=100,10,167.7,5.3,0.05
#记录每次实验样本的置信区间
records=[]
for i in range(g):
    data=norm.rvs(loc=167.7,scale=5.3,size=n)
    #每次样本的值均和标准误
    mean,xsem=data.mean(),sem(data)
    #1-a置信区间CI为：
    records.append([mean-t.isf(a/2,n-1)*xsem,mean+t.isf(a/2,n-1)*xsem])
#置信区间包含住总体μ的概率
count=0
for i in range(g):
    if records[i][0]<μ<records[i][1]:
        count+=1
print('概率为：',count/g)

概率为：0.96

xl,xu=[i[0] for i in records],[i[1] for i in records]
xticks=[i for i in range(len(records))]
plt.vlines(xticks,xl,xu,colors='b',alpha=0.8)
plt.hlines(μ,0,100,'r')

假设检验

t检验

注意！Scipy中的t检验均是双尾检验，若要单尾检验，仅需将p-value/2

单样本t检验

#n=36,xmean=130.83,s=25.74,μ0=140
#H0:μ=μ0=140，即从事铅作业男性工人与正常成年男性的血红蛋白含量均数相等
#H1:μ！=μ0=140，即从事铅作业男性工人与正常成年男性的血红蛋白含量均数不等
#待测样本来自正态总体
ttest_1samp(norm.rvs(size=36,loc=130.83,scale=25.74),140)

Ttest_1sampResult(statistic=-2.6308725605833954, pvalue=0.012577747566433572)

配对样本t检验

#配对实验设计（关联样本）
#差值服从正态总体
#H0:μd=0，即两种方法的测定结果相同
#H1:μd！=0，即两种方法的测定结果不同
data=np.array([[0.84,0.591,0.674,0.632,0.687,0.978,0.75,0.73,1.20,0.87],
               [0.58,0.509,0.5,0.316,0.337,0.517,0.454,0.512,0.997,0.506]])
data

#方法1：相关样本t检验
ttest_rel(data[0],data[1])

Ttest_relResult(statistic=7.925975721054968, pvalue=2.3839523690926263e-05)

#方法2：转化为总体均值为0的单样本t检验
#因p<=2.384,按a=0.05水准，拒绝H0，接受H1，差异有统计学意义。结合题目可以认为这两种方法对脂肪含量的测定结果不同，其中哥特里-罗紫法测定结果较高。
ttest_1samp(data[0]-data[1],0)

Ttest_1sampResult(statistic=7.925975721054968, pvalue=2.3839523690926263e-05)

独立样本t检验

#完全随机设计实验
#两样本均来自正态总体
#两样本的方差相等
#H0:μ1=μ2，即阿卡波糖胶囊组与拜唐苹胶囊组空腹血糖下降值的总体均数相等
#H1:μ1！=μ2，即阿卡波糖胶囊组与拜唐苹胶囊组空腹血糖下降值的总体均数不相等
#待测两样本均来自正态总体
data=np.array([[-0.7,-5.6,2,2.8,0.7,3.5,4,5.8,7.1,-0.5,2.5,-1.6,1.7,3,0.4,4.5,4.6,2.5,6,-1.4],
             [3.7,6.5,5,5.2,0.8,0.2,0.6,3.4,6.6,-1.1,6,3.8,2,1.6,2,2.2,1.2,3.1,1.7,-2]])
data.shape

(2, 20)

#equal_var参数表示，两样本总体方差是否相等，默认为True
#方法1：现有两组样本，直接计算
ttest_ind(data[0],data[1])

Ttest_indResult(statistic=-0.641871949031214, pvalue=0.5248097059442389)

#方法2：现有两组样本的均值和标准差
ttest_ind_from_stats(data[0].mean(),data[0].std(),19,data[1].mean(),data[1].std(),19,equal_var=True)
#p=0.53，按a=0.05水准，不拒绝H0，差异无统计学意义，尚不能认为阿卡波糖胶囊与拜唐苹胶囊对空腹血糖的降糖效果不同。

Ttest_indResult(statistic=-0.641871949031214, pvalue=0.5250226836502808)

正态检验

偏度skewtest检验

skewtest([100, 100, 100, 100, 100, 100, 100, 101])
SkewtestResult(statistic=3.5717766638478086, pvalue=0.0003545677202816322)
skewtest([2, 8, 0, 4, 1, 9, 9, 0])

SkewtestResult(statistic=0.44626385374196975, pvalue=0.6554066631275459)

峰度kurtosistest检验

kurtosistest(list(range(20)))
KurtosistestResult(statistic=-1.7058104152122062, pvalue=0.08804338332528348)
s = np.random.normal(0, 1, 1000)
kurtosistest(s)

KurtosistestResult(statistic=0.35408712426705563, pvalue=0.7232735911597703)

正态分布normaltest检验

a = np.random.normal(0, 1, size=1000)
b = np.random.normal(2, 1, size=1000)
x = np.concatenate((a, b))
x.shape
(2000,)
k2,p=normaltest(x)
#k2为偏度平方＋峰度平方，p为卡方概率
#其中p<0.001,按检验水平a=0.001，拒绝H0，差异有统计学意义，可以认为该样本与正态分布不同。
print('p={:}'.format(p))

p=4.814764728135897e-14

QQ图

#可以看出当df足够大时，t分布近似正态分布
probplot(t.rvs(df=49,size=100),dist='norm',fit=True,plot=plt)

两样本F方差检验

两样本均来自总体总体

#F=var1/Var2 df1=n1-1,df2=n2-1
var1,var2=3.0601**2,2.4205**2
df1,df2=20-1,20-1
F=var1/var2
F

1.5983101734517282

f.sf(F,df1,df2)

0.157643797395044

#因为1.60<2.168,故p>0.1水准，不拒绝H0，差异无统计学意义。尚不能认为阿卡波糖胶囊与拜唐苹胶囊组空腹血糖下降值的总体方差不等。
f.isf(0.05,df1,df2)

2.168251601406261

两样本不服从正态分布

#执行Levene检验以得出均等的方差
a = [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]
b = [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]
c = [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]
stat,p=levene(a,b,c)
#较小的p值表明总体没有相同的方差
p

0.002431505967249681