Python 数学建模——假设检验

前言

  假设检验是概率论中相当重要的内容。一般是先提出一个原假设 $H_0$ 和一个对立的备择假设 $H_1$，通过数学方式获取接受原假设的概率，从而判断是否接受原假设。
  数学建模中有时也需要使用到假设检验，这里主要介绍一下关于正态总体均值的假设检验方法，以及无需正态分布的 KS 检验。
  假设检验大体分为参数假设检验和非参数假设检验两种类型

• 参数检验：总体的分布函数形式已知，但其中含有未知参数。提出关于参数的假设，然后基于样本判定是否接受这个假设。
• 非参数检验：总体的分布函数完全未知。

参数假设检验

单个总体均值的假设检验

用途：已知一个量 $X$ 服从正态分布，我现在观测到了一些 $X$ 的值，那么能不能认为 $X$ 的均值是 ${\mu }_0$。

  随机变量 $X\sim N(\mu ,{\sigma }^2)$，并且其中的 $\mu$ 是未知的，是我们要假设检验的参数。并且已经有 $n$ 个随机变量的观测值 $X_1,\cdots ,X_n$。提出假设：

• 原假设 $H_0:\mu ={\mu }_0$。
• 备择假设 $H_1:\mu \ne {\mu }_0$。

  下面基于两种情况，采用不同的方式进行假设检验：

$\sigma$ 已知

  检验统计量 $Z=\frac{\overline{X}-{\mu }_0}{\sigma /\sqrt{n}}$，检验的显著性水平为 $\alpha$，标准正态分布的上 $\alpha /2$ 分位数记为 $z_{\alpha /2}$。
  若 $|Z|>z_{\alpha /2}$，那么拒绝原假设，接受备择假设；否则接受原假设。

$\sigma$ 未知

  检验统计量 $T=\frac{\overline{X}-{\mu }_0}{s/\sqrt{n}}$。这个 $S$ 是根据 $n$ 个样本算出来的标准差。检验的显著性水平为 $\alpha$，自由度为 $n-1$ 的 $t$ 分布的上 $\alpha /2$ 分位数记为 $t_{\alpha /2}(n-1)$。
  若 $|T|\ge t_{\alpha /2}(n-1)$，那么拒绝原假设，接受备择假设；否则接受原假设。

一般 $\sigma$ 是未知的，即使已知也能当做未知。下面是 $\sigma$ 未知的代码：

from scipy.stats import ttest_1samp

# 在这里输入要检验的总体 data
data = [1,1,4,5,1,4,1,9,1,9,8,1,0]

# 即我们期望的平均水平 mu_0
mu = 3

tset, pval = ttest_1samp(data, mu)
print(pval)

  给出来的pval也称作 $p$ 值，一般 $p>0.05$ 即可接受原假设。

两个总体均值的假设检验

用途：已知两个量 $X$，$Y$ 服从正态分布，我现在观测到了一些 $X$ 的值和 $Y$ 的值，那么能不能认为 $X$ 和 $Y$ 的均值差不多。（$X$，$Y$ 需要方差相同，可以放在论文假设里面）

  两个总体均值的假设检验的前提是两个正态分布的方差均已知且相等，即 $X\sim N({\mu }_1,{\sigma }^2)$，$Y\sim N({\mu }_2,{\sigma }^2)$，并且两个随机变量相互独立。已经获取了 $X$ 的 $n_1$ 个观测值和 $Y$ 的 $n_2$ 个观测值，提出假设：

• 原假设 $H_0:{\mu }_1={\mu }_2$。
• 备择假设 $H_1:{\mu }_1\ne {\mu }_2$。

  检验的显著性水平为 $\alpha$，检验统计量：
$$T=\frac{\overline{X}-\overline{Y}}{\sqrt{\frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-1}}\cdot \sqrt{\frac{1}{n_1}+\frac{1}{n_2}}}$$

若 $|T|\ge t_{\alpha /2}(n_1+n_2-2)$，那么拒绝原假设，接受备择假设；否则接受原假设。

接受原假设代表两个随机变量的分布，在显著性水平为 $\alpha$ 的情况下接受它们没有显著的差异。
启发：如果题目问两个量是否有显著差异，刚好我们又可以认为两个量均服从正态分布，方差相等并且它们相互独立，那么可以采用假设检验的方法验证。

参考代码

import numpy as np
from statsmodels.stats.weightstats import ttest_ind
from scipy.stats import levene

# a，b 数据量可以不同
a=np.array([1,1,4,5,1,4])
b=np.array([1,9,1,9,8,1,0])
tstat, pvalue, df=ttest_ind(a, b, value=0,usevar='pooled' if (levene(a,b).pvalue>0.05) else 'unequal')
print(' 检验统计量为：', tstat)
print('p 值为：', pvalue)
print(' 自由度为：', df)

非参数假设检验

  在实际建模中，对样本数据服从什么分布完全是未知的，需要进行非参数假设检验。下面是两种非参数假设检验方法。

分布拟合检验——卡方检验

用途：很少用到，特别是这个 $k$ 取几，区间的边界是多少，感觉很难说清楚。

  假设给了随机变量的 $n$ 个观测值。那么 ${\chi }^2$ 检验的步骤是：

1. 建立待检验假设 $H_0$：样本总体对应的随机变量 $X$ 的分布函数是 $F(x)$。
2. 将 $\mathbb{R}$ 分成 $k$ 个区间，令 $p_i$ 为 $X$ 落入第 $i$ 个区间的概率，$m_i$ 是给的样本中落入第 $i$ 个区间的频数。

这个 $k$ 没说具体要怎么取，感觉可以自己控制一下，哪个效果好就用哪个。
按照这个逻辑，我们直观的感受就是，$m_i$ 越接近 $n{p}_i$，说明假设给出的分布函数 $F(x)$ 效果就越好。

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)$ 中未知参数的个数。
4. 对于给定的显著性 $\alpha$，确定卡方分布的上 $\alpha$ 分位数，即确定 ${\chi }_{\alpha }^2$ 使得 P { χ 2 ( k − 1 − r ) > χ α 2 } = α P\{{{\chi }^{2}}(k-1-r)>{{\chi }_{\alpha }^{2}}\}=\alpha P{χ2(k−1−r)>χα2​}=α。
5. 当 3 中算出来的 ${\chi }^2$ 小于 4 中算出来的 ${\chi }_{\alpha }^2$ 时，接受假设 $H_0$。否则拒绝这个假设。

KS 检验（Kolmogorov-Smirnov 检验）

用途：当你猜测某个变量 $X$ 服从一种特定的分布的时候，可以使用 KS 检验验证一下。

  样本可以得到一个经验分布函数 $F_n(x)$，我们认为样本对应的随机变量服从一个理论分布函数 $F(x)$。用检验量 $T={\sup }_x|F_n(x)-F(x)|$ 来表示两个分布函数之间的差距。
  Python 包根据这些样本使用一些方式（书里没有说清楚）得到一个值 $t$，并获取 $p$ 值 p = P { T ≥ t } p=P\{T\ge t\} p=P{T≥t}，当这个 $p$ 值很小的时候（一般是 $p<0.05$）就拒绝假设，即认为所指定的分布是不可接受的。
  下面是代码示例，其中scipy.stats.kstest可以用来做 KS 检验。

#程序文件Pex4_21.py
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as ss
a=np.array(pd.read_excel("Pdata4_6_1.xlsx"))
w=a[:,1::2]; w=w.flatten()
# 到这里 w 已经变成了一维数组，里面是体重的数据
mu=w.mean(); s=w.std(ddof=1)  #计算样本均值和标准差
print("均值和标准差分别为：", (mu, s))
statVal, pVal=ss.kstest(w,'norm',(mu,s))
print("统计量和P值分别为：", [statVal, pVal])

ketest更多用法参考文献：Python SciPy stats.kstest用法及代码示例 - 纯净天空

Wilcoxon 检验

  分为 Wilcoxon 符号秩检验和 Wilcoxon 秩和检验，两者的原理和用法是不相同的。

Wilcoxon符号秩检验

用途：比较两个配对序列之间的差异。换句话说，也就是判断两个样本的中位数是否有显著差异。

import scipy.stats as ss

# x,y 分别是这两个配对序列
x=[310,350,370,377,389,400,415,425,440,295,325,296,250,340,298,365,375,360,385]
y=[320]*len(x)

# 重要参数：
# correction: 默认 False，小样本需要 True。对于小样本没有明确定义。应该可以认为 < 40 是小样本
# alternative: 从['greater','less','two-sided']中选择，表示备择检验中 x 的中位数与 y 中位数的关系。
#    - 'greater'原假设： x 中位数不显著大于 y 中位数
#    - 'less'原假设： x 中位数不显著小于 y 中位数
#    - 'two-sided'原假设： x 中位数不显著大于也不显著小于 y 中位数
print(ss.wilcoxon(x,y,correction=True,alternative='greater'))

  打印出来的结果中有一个 $p$ 值，$p$ 值越大两个样本的中位数越一致，$p>0.05$ 即可认为没有显著差异。
参考文献：非参数统计的Python实现—— Wilcoxon 符号秩检验_非参数 wilcoxon 符号秩检验-CSDN博客

Wilcoxon秩和检验

用途：用于比较两组数据是否服从同一个分布。不需要配对，并且每组数据都是对同一个属性的不同采样值。

import scipy.stats as ss

# 两组列表数据
data1 = [1,1,4,5,1,4,1,9,1,9]
data2 = [8,1,0,1,1,4,11,4]

stat, p = ss.ranksums(data1, data2)
print(p)

  同样打印出来的结果中有一个 $p$ 值，$p$ 值越大越有把握认为服从同一分布，$p>0.05$ 即可认为服从同一分布。
参考文献：【统计分析】Python实现威尔克逊秩和检验（Wilcoxon Rank Sum test）_wilcoxon秩和检验python-CSDN博客