区间估计、假设检验、经验分布、Q-Q图、非参数假设检验（卡方、柯尔莫哥洛夫、秩和）、Matlab实现案例

区间估计

例1

• 分别使用金球和铂球测定引力常数（单位：$10^{-11}{m}^3\cdot k{g}^{-1}\cdot s^{-2}$）。

（1）用金球测定观察值为6.683，6.681，6.676，6.678，6.679，6.672。
（2）用铂球测定观察值为6.661，6.661，6.667，6.667，6.664。
设测定值总体为$N(\mu ,{\sigma }^2),\mu ,{\sigma }^2$均为未知，试就（1）,（2）两种情况分别求$\mu$的置信度为0.9的置信区间，并求${\sigma }^2$的置信度为0.9的置信区间。
（3）设用金球和铂球测定时总体的方差相等，求两个测定值总体均值差的置信度为0.90的置信区间。
解：
（1）$\mu ,{\sigma }^2$均未知时，$\mu$的置信度为0.9的置信区间为
$$\left(\bar{X}-\frac{S}{\sqrt{n}}{t}_{\alpha /2}(n-1),\bar{X}+\frac{S}{\sqrt{n}}{t}_{\alpha /2}(n-1)\right)$$
（2）$\mu ,{\sigma }^2$均未知时，${\sigma }^2$的置信度为0.9的置信区间为
$$\left(\frac{(n-1)S}{{\chi }_{\alpha /2}^2(n-1)},\frac{(n-1)S}{{\chi }_{1-\alpha /2}^2(n-1)}\right)$$
（3）总体均值差的置信度为0.90的置信区间为
$$\left(\overline{X_1}-\overline{X_2}\pm t_{\alpha /2}(n_1+n_2-2)S_w\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}\right)$$

clc,clear
x1=[6.683,6.681,6.676,6.678,6.679,6.672];
x2=[6.661,6.661,6.667,6.667,6.664];
%如果检验在 5%(默认) 的显著性水平上拒绝原假设，则结果 h 为 1，否则为 0。
%p假设检验的p值。
%均值的置信区间 ci，以及包含检验统计量信息的结构体 stats。
[h1,p1,ci1,st1]=ttest(x1,mean(x1),'Alpha',0.1)%t检验可求得均值的置信区间
[h2,p2,ci2,st2]=ttest(x2,mean(x2),'Alpha',0.1)
[h3,p3,ci3,st3]=vartest(x1,var(x1),'Alpha',0.1)
[h4,p4,ci4,st4]=vartest(x2,var(x2),'Alpha',0.1)
[h,p,ci,st]=ttest2(x1,x2,'Alpha',0.1)%双样本t检验

# 方差既可以使用numpy函数，也可以使用pandas函数。
# 
# numpy 中计算的方差就是样本方差本身：
# 
#         使用场景为：拥有所有数据的情况下，计算所有数据的标准差时使用，即最终除以n，而非n-1
# 
# 
# pandas 中计算的方差为无偏样本方差：
# 
#         使用场景为：只有部分数据但需要求得总体的标准差时使用，当只有部分数据时，根据统计规律，除以n时计算的标准差往往偏小，因此需要除以n-1，即n-ddof。
# 
# 由于是用于样本数据，所以采用pandas的方差函数。
import pandas as pd
import numpy as np
from scipy import stats
x1 = pd.Series([6.683,6.681,6.676,6.678,6.679,6.672])
x2 = pd.Series([6.661,6.661,6.667,6.667,6.664])
x1_mean,x1_var,x1_std=stats.bayes_mvs(x1,alpha=0.90)
x2_mean,x2_var,x2_std=stats.bayes_mvs(x2,alpha=0.90)
print('总体x1均值的置信度为0.90的置信区间是{}'.format(x1_mean))
print('总体x1方差的置信度为0.90的置信区间是{}'.format(x1_var))
print('总体x2均值的置信度为0.90的置信区间是{}'.format(x2_mean))
print('总体x2方差的置信度为0.90的置信区间是{}'.format(x2_var))

注：python区间估计和假设检验部分可查看我的这份笔记。
语雀：https://www.yuque.com/chenyujiao-4zrlp/sxmd8u/adwx4l
jupter notebook：http://localhost:8888/notebooks/Documents/python%E6%95%B0%E6%8D%AE%E5%88%86%E6%9E%90/%E7%AC%AC6%E7%AB%A0%20%E7%AE%80%E5%8D%95%E7%BB%9F%E8%AE%A1%E6%8E%A8%E6%96%AD.ipynb

经验分布函数

设$x_1,x_2,\dots ,x_n$是总体$F$的一个容量为$n$的样本值。从小到大排序后：$x_{(1)}\le x_{(2)}\le \dots \le x_{(n)}$
$$F_n(x)=\{\begin{array}{ll}0,& x<x_{(1)}\\ \frac{k}{n},& x_{(k)\le x<x_{(k+1)}}\\ 1,& x\ge x_{(n)}\end{array}$$
对于经验分布函数，格里汶科（Glivenko）证明了当$n\to +\infty$时，$F_n(x)$以概率1一致收敛于总体分布函数$F(x)$。

例2

下面给出了84个伊特拉斯坎任男子的头颅的最大宽度（mm），计算经验分布函数并画出经验分布函数图形。
141 148 132 138 154 142 150 146 155 158 150 140 147 148 144 150 149 145 149 158 143 141 144 144 126 140 144 142 141 140 145 135 147 146 141 136 140 146 142 137 148 154 137 139 143 140 131 143 141 149 148 135 148 152 143 144 141 143 147 146 150 132 142 142 143 153 149 146 149 138 142 149 142 137 134 144 146 147 140 142 140 137 152 145

clc,clear
load 'data.mat' #数据在一维数组data中
[ycdf,xcdf,n]=cdfcalc(data);%注意得出的经验分布的值没有保留重复值
cdfplot(data)
title("")
hold on
plot(xcdf,ycdf(2:end),'r*')

"""经验分布函数"""
import numpy as np
from statsmodels.distributions.empirical_distribution import ECDF
import pandas as pd
data=pd.read_excel(r'ex7_5.xls')
ecdf=ECDF(data['data'])
sort_data=data.sort_values(by='data')['data']
ecdf(sort_data)

import matplotlib.pyplot as plt
plt.figure()
plt.step(sort_data,ecdf(sort_data))
plt.plot(sort_data,ecdf(sort_data),'r-')

Q-Q图

作用

• Q-Q图(Quantile-Quantile Plot)是检验拟合优度的好办法。对于一组观察数据$x_1,x_2,...,x_n$，利用参数估计方法确定了分布模型的参数$\theta$后，分布函数$F(x;\theta )$就知道了，现在希望知道的是观测数据与分布模型的拟合效果如何。

原理

• Q-Q图的基本思想是将经验分布函数的分位数点和分布模型的理论分布分位数点作为一堆数组画在直角坐标图上，就是一个点，$n$个观测数据对应$n$个点，如果这$n$个点看起来像一条直线，说明观测数据与分布模型的拟合效果很好。

计算步骤

1. 将$x_1,x_2,...,x_n$依大小顺序排列成$x_{(1)}\le x_{(2)}\le \cdots \le x_{(n)}$。
2. 取$y_i=F^{-1}((i-1/2)/n),i=1,2,\dots ,n$。
3. 将$(y_i,x_{(i)}),i=1,2,\dots ,n$，这$n$个点画在直角坐标图上。
4. 如果这$n$个点看起来呈一条$45°$角的直线，则拟合效果很好。

续例2

• 如果这些数据来自于正态总体，求该正态分布的参数，试画出它们的Q-Q图，判断拟合效果。

clc,clear
load 'data.mat'
mean_data=mean(data);
std_data=std(data);
%PD = ProbDistUnivParam（DistName， Params） 创建 PD，这是一个 ProbDistUnivParam 对象，它表示概率分布。此分布由 DistName 指定的参数分布定义，参数由数值矢量参数指定。
pd=makedist('normal','mu',mean_data,'sigma',std_data);
%qqplot（X） %检查数据X是否副总正态分布，这和正态概率函数差不多
%qqplot（X，Y） %检验数据X与Y是否服从相同的分布（他们的每个分位数是否都相近相似）
qqplot(data,pd)

from scipy import stats
data.mean()
data.std()
stats.probplot(sort_data, dist="norm",plot=plt)
plt.show()

非参数检验

${\chi }^2$拟合优度检验

步骤

1. 建立待检验假设$H_0$:总体$X$的分布函数为$F(x)$。
   1. 若总体$X$是离散型的，建立待检假设$H_0$:总体$X$的分布律为 P { X = x i } = p i , i = 1 , 2 , . . . P\{X=x_i\}=p_i,i=1,2,... P{X=xi​}=pi​,i=1,2,...。
   2. 若总体$X$是连续型的，建立待检假设$H_0$:总体$X$的概率密度为$f(x)$。
2. 将数轴分成$k$个区间，令$p_i$为分布函数$F(x)$的总体$X$在第$i$个区间内取值的概率，设$m_i$为$n$个样本观察值落入第$i$个区间上的个数，也称为组频数。
3. 选取统计量${\chi }^2={\sum }_{i=1}^k\frac{(m_i-n{p}_i{)}^2}{n{p}_i}={\sum }_{i=1}^k\frac{m_i^2}{n{p}_i}-n$，如果$H_0$为真，则${\chi }^2\sim {\chi }^2(k-1-r)$，其中$r$为分布函数$F(x)$中未知参数的个数，$k$为分组的组数。
4. 给定显著性水平$\alpha$，确定${\chi }_{\alpha }^2$并计算样本统计量${\chi }^2$的观察值。
5. 作出判断：若${\chi }^2<{\chi }_{\alpha }^2$，则接受$H_0$；否则拒绝$H_0$，即不能认为总体$X$的分布函数为$F(x)$。

续例2

• 试检验这些数据是否来自正态总体（取$\alpha =0.1$）。

解：$H_0$：头颅的最大宽度$X$服从正态分布$N(143.7738,5.9705^2)$。

${\chi }^2=88.3408-88=4.3408<{\chi }_{\alpha }^2(k-r-1)={\chi }_{0.1}^2(7-2-1)=7.7794$
所以接受假设，认为这些数据是来自正态总体。

clc,clear
load 'data.mat'
%创建正态分布的方法一：这种方法两个参数是待估的
pd_f=fitdist(data','Normal')%对数据进行概率分布对象拟合，用数据拟合正态分布得到的参数的点估计和区间估计
%创建正态分布的方法二：这种方法两个参数是确定的
mean_data=mean(data);
std_data=std(data);
pd_m=makedist('normal','mu',mean_data,'sigma',std_data);%创建概率分布对象，有和fitdist一样的方法

%得到分布函数的方法一：
pd_cdf=@(x) cdf(pd_f,x)或pd_cdf=@(x) cdf(pd_m,x)
%得到分布函数的方法二：
pd_cdf=@(x) normcdf(x,mean_data,std_data)
%qqplot（X） %检查数据X是否服从正态分布，这和正态概率函数差不多
%qqplot（X，Y） %检验数据X与Y是否服从相同的分布（他们的每个分位数是否都相近相似）
qqplot(data,pd)

%默认bins为10，会使得观测数均在5个以上。
%‘emin’非负整数，默认值5，指定一个区间对应的最小理论频数
%‘frequency’与x等长的向量， 指定x中各元素出现的频数
%‘alpha’默认值0.05 ，指定检验的显著性水平
[h,p,st]=chi2gof(data,'cdf',pd_cdf,'NParams',2,'Alpha',0.1) %cdf指分布函数；‘expected’指定各区间的理论频数，与‘cdf’不能同时出现

[h,p,st]=chi2gof(data,'cdf',pd_f,'Alpha',0.1) %由fitdist得到的分布对象，r=2,统计量的自由度为k-r-1=4
[h,p,st]=chi2gof(data,'cdf',pd_m,'Alpha',0.1) %由makedist得到的分布对象，r=0,统计量的自由度为k-r-1=6

%本题中，两个参数是估计出来的，所以选第一种或第二种，r=2才对

结果：

h =0
p =0.3618
st = 

包含以下字段的 struct:

chi2stat: 4.3408 %样本统计量的结果
df: 4 %自由度
edges: [126.0000 135.6000 138.8000 142.0000 145.2000 148.4000 151.6000 158.0000] %分组的边界值
O: [7 7 22 15 15 10 8]%实际频数
E: [7.1816 9.8204 15.1865 17.7409 15.6563 10.4375 7.9768]%理论频数

>>> chi2inv(0.9,st.df)%临界值
ans =
    7.7794

例3

• 在一批灯泡中抽取300只作寿命试验，取$\alpha =0.05$，试检验假设：$H_0$：灯泡寿命服从指数分布。

$$f(t)=\{\begin{array}{ll}0.005e^{-0.005t},& t\ge 0\\ 0,& t<0\end{array}$$

寿命$t/h$	$0<t\le 100$	$100<t\le 200$	$200<t\le 300$	$t>300$
灯泡数	121	78	43	58

解：

clc,clear
edges=[0:100:300 inf];
bins=[50 150 250 inf];
num=[121 78 43 58];
pd=makedist('exp',200);
pd_cdf=@(x) cdf(pd,x);
[h,p,st]=chi2gof(bins,'Edges',edges,'cdf',pd_cdf,'Frequency',num,'NParams',0)
chi2=chi2inv(0.95,st.df)

结果：
h =     0
p =    0.6052
st = 
  包含以下字段的 struct:

    chi2stat: 1.8450
          df: 3
       edges: [0 100 200 300 Inf]
           O: [121 78 43 58]
           E: [118.0408 71.5954 43.4248 66.9390]
chi2=		7.8147

从h=0和chi2stat<chi2和p>0.05都可以说明接受原假设，即认为这批灯泡寿命服从参数为200的指数分布。

柯尔莫哥洛夫检验

原理

• $\chi 2$拟合优度检验实际上是检验$p_i=F_0(a_i)-F(a_{i-1})=p_{i0}(i=1,2,...,k)$的正确性，并未直接检验原假设的分布函数$F_0(x)$的正确性，柯尔莫哥洛夫检验直接针对原假设$H_0:F(x)=F_0(x)$，这里分布函数$F_{(}x)$必须是连续型分布。柯尔莫哥洛夫检验基于经验分布函数作为检验统计量，检验理论分布函数与样本分布函数的拟合优度。
• 设总体$X$服从连续分布，$X_1,X_2,...,X_n$是来自总体$X$的简单随机样本，$F_n(x)$为经验分布函数。定义

$$D_n=\underset{-\infty <x<+\infty }{\sup }|F_n(x)-F(x)|$$
根据大数定律，当$n$区域无穷大时，$D_n$依概率收敛到0.

步骤

1. 原假设和备择假设

$$H_0:F(x)=F_0(x),H_1:F(x)\ne F_0(x)$$

2. 选取检验统计量

$$D_n=\underset{-\infty <x<+\infty }{\sup }|F_n(x)-F(x)|$$
当$H_0$为真时，$D_n$有偏小趋势，则拟合得越好。

3. 确定拒绝域。给定显著性水平$\alpha$，查$D_n$极限分布表，求出$t_{\alpha }$满足 P { n D n ≥ t α } = α P\{\sqrt n D_n\ge t_{\alpha}\}=\alpha P{n ​Dn​≥tα​}=α，即拒绝域为$[t_{\alpha }/\sqrt{n},+\infty ]$。

$\alpha$	0.1	0.05	0.01
$t_{\alpha }$	1.22	1.36	1.63

4. 作判断。计算统计量得观察值，如果$D_n$落在拒绝域则拒绝原假设。

续例2

• 试用柯尔莫哥洛夫检验法检验这些数据是否服从正态分布($\alpha =0.05$)。

解：

1. $H_0:X\sim N(143.7738,5.9705^2)$
2. 拒绝域是 { t 0.05 / n , + ∞ } = { 0.1484 , + ∞ } \{t_{0.05}/\sqrt n,+\infty\}=\{0.1484,+\infty\} {t0.05​/n ​,+∞}={0.1484,+∞}
3. 计算经验分布函数$F_n(x_i)$和理论分布函数值$F(x_i)$，并计算统计量

$$D_n={\sup }_{x_i}|F_n(x)-F(x)|=0.0851<0.1484$
所以接受原假设，即认为服从正态分布。

clc,clear
load 'data.mat'
%% 方法一：由该方法的步骤计算
mean_data=mean(data);
std_data=std(data);
[ycdf,xcdf,n]=cdfcalc(data);%经验分布
ycdf=ycdf(1:end-1);%% 这里有一点疑惑，按照经验分布函数的定义，不应该是去掉第一个数据0吗
pd_cdf=@(x) normcdf(x,mean_data,std_data);%理论分布
%pd_cdf=normcdf(xcdf,mean_data,std_data);
Dn=max(abs(ycdf - pd_cdf(xcdf)));
if Dn*sqrt(length(data))<1.36
    disp("接受原假设")
else
    disp("拒绝原假设")
end

%% 方法二：调用工具箱
pd=makedist('Normal','mu',mean_data,'sigma',std_data);
[h,p,st,cv]=kstest(data,'CDF',pd,'Alpha',0.05)%cv是临界值，和上面结果不一样是因为这里的t更准确

结果：

h =    0 %代表接受原假设
p =    0.5484 %p值大于显著性水平0.05，所以接受原假设
st =    0.0851 %代入样本数据得到的统计量的值，大于临界值cv，所以接受原假设
cv =    0.1461 %临界值

秩和检验

作用

• 检验假设$H_0$：两个总体$X$与$Y$有相同的分布。

步骤

设分别从总体$X、Y$两个总体独立抽取大小为$n_1、n_2$，且$n_1\le n_2$。

1. 将两个样本混合起来，按照数值从小到大统一排序，每个数据对应的序数为秩。
2. 计算取自总体$X$的样本所对应的秩之和，用$T$表示。
3. 根据$n_1、n_2$与水平$\alpha$，查秩和检验表，得秩和下限$T_1$和秩和上限$T_2$。
4. 若$T\le T_1或T\ge T_2$，则否定原假设，认为两个总体分布有显著差异。

原理

如果两总体分布无显著差异，那么$T$不应太大或太小。

例3

• 试讨论烘干温度对抗弯强度在水平$\alpha =0.05$下是否有显著影响？

解：

1. 排序

2. $n_1=6,n_2=9,T=39$
3. 查表$T_1=33,T_2=63$
4. $T_1\le T\le T_2$，接受原假设，认为没有显著差异

clc, clear
x=[41.5	42.0	40.0	42.5	42.0	42.2	42.7	42.1	41.4];
y=[41.2	41.8	42.4	41.6	41.7	41.3];
yx=[y,x];
yxr=tiedrank(yx);%计算秩
yr=sum(yxr(1:length(y)));%计算y的秩和

%% 利用工具箱
[p,h,st]=ranksum(y,x,'Alpha',0.05)