Task01|GriModel统计分析（下）|方法论与一元数值检验|假设检验1

假设检验

• 假设检验可以用于验证模型的有效性
  例子：
  • 回归分析中模型参数的显著性检验（基于统计模型）
  • 时间序列分析中单位根检验、白噪声检验

假设检验两大类（模型）

• 基于统计模型的假设检验
• 不基于统计模型的假设检验
  直接从数据出发，检验数据的的某些性质
  exp:正态性检验，双样本t检验，方差检验

假设检验的原理

根据样本信息与已知信息，对一个描述总体性质的命题进行数理上的检验并分析

假设检验两大类||(参数)

• 参数假设检验
  假设是关于总体的一个参数或者参数的集合
  exp:总体均值
• 非参数假设检验
  假设无法用一个参数集合表示
  exp:正态性检验

假设检验的推导

1. 假设的建立
   exp：样本服从正态分布，抽取样本进行参数假设检验之样本总体均值检验.
   将建设检验的的命题化成两个互斥假设，原假设H0和备择假设H1（通过判断是否拒绝假设来判断是否接受假设）
   通常将想要验证的命题作为备择假设H1
2. 单参数的三种类型与原假设的注意事项

• 在检验某个样本所在的通体是否服从某个特定的分布的时候这样设置:
  H0:样本所在的总体服从某分布
  H1:样本所在的总体不服从某分布

• 在单参数检验中‘=’只会出现在原假设H0中，而不会出现在备择假设H1中

三种参数检验问题

• 均值检验:
  $$\begin{gathered} H_0:\mu ⩽{\mu }_0\leftrightarrow H_1:\mu >{\mu }_0（单边检验） \\ {H}_0:\mu ⩾{\mu }_0\leftrightarrow H_1:\mu <{\mu }_0（单边检验） \\ {H}_0:\mu ={\mu }_0\leftrightarrow H_1:\mu \ne {\mu }_0（双边检验） \end{gathered}$$

另一种表达方式：

$$\begin{gathered} H_0:\mu ={\mu }_0\leftrightarrow H_1:\mu >{\mu }_0 \\ {H}_0:\mu ={\mu }_0\leftrightarrow H_1:\mu <{\mu }_0 \\ {H}_0:\mu ={\mu }_0\leftrightarrow H_1:\mu \ne {\mu }_0 \end{gathered}$$

在以后的参数检验中，我们统一将原假设的符号设置为=

临界值，拒绝域，显著性水平

Exp1
样本均值$\bar{x}$
总体均值$\mu$
有参数假设检验
$H_0:\mu ⩽110\leftrightarrow H_1:\mu >110$
转换为:
$H_0:\mu \ne 110\leftrightarrow H_1:\mu =110$

样本均值$\bar{x}$是总体均值$\mu$的无偏估计
如果$\mu >110$那么$\bar{x}$很有可能大于110
如果实际样本$\bar{x}\gg 110$那么原假设很有可能不成立

临界值：
假设一个临界值C,如果满足$\bar{x}-110>C$则拒绝原假设

拒绝域
{ x ˉ    : x ˉ > 110 + C } \left\{ \bar{x}\,\,: \bar{x}>110+C \right\} {xˉ:xˉ>110+C}
样本计算结果落入拒绝域中拒绝原假设，不同假设检验的拒绝域不尽相同，核心的逻辑是完全相同的

• 如何确定临界值C?
  概率来确定
  在检验中有两种犯错类型:

1. 拒真(A)
   拒真概率$\alpha$
2. 取伪(B)
   取伪概率$\beta$

两种犯错概率相互对立，犯了错要么犯A要么犯B

在样本量给定的条件下，如果降低某一类错误必然导致另一类错误概率的上升

调整假设检验规则可以降低某一类错误概率

通常的做法是只限制犯第一类错误(拒真)的概率

显著性水平$\alpha$

临界值C要保证发生第一类错误的概率$\alpha$需要在一个较小的范围($\alpha =0.05/0.1$)

确定临界值C的标准是：原假设$H_0$为真,数据落入拒绝域的概率恰好为$\alpha$(拒真的概率为$\alpha$)

在Exp1中，这个概率可以写为:

$P_{H_{0}istrue}\left(\bar{x}-{\mu }_0>C\right)=P\left(\bar{x}-110>C\right)=\alpha$

用分布处理概率——构造检验统计量

$P\left(\bar{x}-110>C\right)$
将上式子构造成某个分布的分位数形式

• 此处的110是在原假设$H_0$正确下的值
  $E(\bar{x})=\mu$
  概率实际为:
  $P\left(\bar{x}-E\left(\bar{x}\right)>C\right),E\left(\bar{x}\right)={\mu }_0=110$

• $\bar{x}$服从正态分布

可以构造t统计量
$P(\bar{x}-E(\bar{x})>C)=P(\frac{\bar{x}-E(\bar{x})}{s}>\frac{C}{s})=\alpha ,\sim \frac{\bar{x}-E(\bar{x})}{s}$

这意味着，$\frac{C}{s}$恰好是$t_{n-1}(1-\alpha )$分位点，

对于一个确定的分布而言,分位点是已知的

可求得分位点:
$C=s\cdot t_{n-1}(1-\alpha )$

我们将C分位点代入上述式子中有
$P(\frac{\bar{x}-E(\bar{x})}{s}>\frac{C}{s})=P(\frac{\bar{x}-{\mu }_0}{s}>\frac{s\cdot t_{n-1}(1-\alpha )}{s})=P(\bar{x}>{\mu }_0+s\cdot t_{n-1}(1-\alpha ))=\alpha$

只要$\bar{x}>110+s\cdot t_{n-1}(1-\alpha )$
可以在显著水平$\alpha$下拒绝原假设

在上面的过程中，我们用已知的统计量构造了一个服从某个特定分布的统计量，构造出来的统计量是检验统计量

假设检验的思路:

1. 构造检验统计量
2. 输出对应分布的分位点
3. 计算临界值
4. 做出判断

检验统计量:
$$t=\frac{\bar{x}-{\mu }_0}{s}$$

对应t分布的分位点为:
$$t_{n-1}(1-\alpha )$$

拒绝域为:
$$\bar{x}>110+s\cdot t_{n-1}(1-\alpha )$$

t分布假设检验

单边检验

import numpy as np
import pandas as pd
from scipy.stats import t
n=25
x_bar=108.2
s=4
mu=110

#计算检验统计量 t值
tvalue=(x_bar-mu)
print('t值为:{}'.format(tvalue))

#输出分位点*******
‘’‘
ppf:单侧左分位点
isf：单侧右分位点
interval:双侧分位点
’‘’
T_isf=t.isf(0.05,n-1)
# 0.05显著性水平，n-1 自由度
#备择假设是大于号 选用单侧右分位点isf
#小于号选择单侧左分位点ppf

print('分位点为:{}'.format(T_isf))

#拒绝域
Deny_domain=110+s*T_isf
print('拒绝域的临界点为:{}'.format(Deny_domain))

##判断
print('样本均值是否位于拒绝域:{}'.format(x_bar>Deny_domain))

#直接用检验统计量与分布分位点判断
# 直接用检验统计量与分布分位点判断
print('检验统计量是否位于拒绝域：{}'.format(tvalue>T_isf))

拒绝域可以直接用检验统计量与对应的分布分位点表示

$t>t_{n-1}(1-\alpha )$

双边检验

两个假设:
$H_0:\mu =110\leftrightarrow H_1:\mu \ne 110$
拒绝域:
$\left|t\right|>\left|t_{n-1}\left(1-\frac{\alpha }{2}\right)\right|$

T_int=t.interval(1-0.0.5,n-1)
#分位点1-alpha,n-1自由度

#统计量t的绝对值
np.abs(tvalue)
（双边检验的重点）

#双侧分位点
T_int

p值（一个不变的判断是否拒绝原假设的指标）

采用拒绝域进行假设检验有个缺点，分位点值与显著性水平$\alpha$是相关的

如果要在不同显著水平检验检验，需要计算不同的分位点再比较

p值是在确定的样本观测值下，给出的能拒绝原假设的最小显著性水平

• p值越小越可以拒绝原假设，p值必须小于置信水平才能拒绝
  如果p为0.001，认为在0.01可以拒绝原假设
  如果p为0.025比0.01的置信水平要大，但小于0.05，认为在0.05置信水平下可以拒绝原假设，在0.01置信水平下不可以拒绝

求p值

· 若$H_1$的符号为$\ne$，则：$pvalue=P\left(\left|X\right|>\left|Teststatistics\right|\right)$

· 若$H_1$的符号为>，则：$pvalue=P\left(X>Teststatistics\right)$

· 若$H_1$的符号为<，则：$pvalue=P\left(X<Teststatistics\right)$

X为服从某个特定分布的变量，p值本质是累积概率

'''
sf:右尾累积概率
cdf:左尾累积概率
'''
# 若备择假设为mu>110
pvalue=t.sf(tvalue,n-1) 
print('备择假设为mu>110的p值为：{}'.format(pvalue))

# 若备择假设为mu<110
pvalue=t.cdf(tvalue,n-1)
print('备择假设为mu<110的p值为：{}'.format(pvalue))

# 若备择假设为mu不等于110
#双边检验
pvalue=t.cdf(tvalue,n-1)*2 # 之所以是左尾累积概率的两倍，是因为右尾累积概率大于0.5，而p值不可能大于1。
print('备择假设为mu不等于110的p值为：{}'.format(pvalue))

#需要筛选掉不合理的值

在后续的学习当中，我们统一使用p值进行假设检验。

1.1.3 假设检验的基本步骤——基于p值⭐⭐※

总结一下基于p值的假设检验基本步骤

1. 确定备择假设$H_1$，备择假设的符号决定了我们使用何种累积概率。
2. 明确检验统计量Test statistics的公式，不同假设检验都有其确定的检验统计量，查阅资料找到它！
3. 明确检验统计量所服从的分布，这样我们才能计算累积概率。
4. 根据备择假设$H_1$与检验统计量Test statistics计算p值。
5. 将p值与显著性水平$\alpha$比较，若$p>\alpha$，则不能拒绝原假设；若$p<\alpha$，则可以拒绝原假设。

1.2假设检验的分类

import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
from IPython.display import display

数值型数据

• 一元数值类型
• 多元数值类型

一元数值类型的假设检验

1. 正态性检验
2. 比较一组数据的总体均值与一个固定值是否相等的检验
3. 比较两组数据的总体均值之间是否相等的检验
4. 比较两组以上的多组数据的总体均值之间是否相等的检验

如果数据大致符合正态分布，那么可以使用参数检验-t检验
如果数据不是正态分布的，那么可以使用一些非参数检验。

2.1正态性检验

参数检验比非参数检验更灵敏
三种判断数据正态性的方法:

1. 可视化判断-正态分布概率图
2. Shapiro-Wilk检验
3. D’Agostino’s K-squared检验

2.1.1概率图

1. 排序
2. 计算排序后数据𝑥对应的分布分位数
3. 分位数为横轴，有序样本值为纵轴
   if 服从目标分布 数据点会近似地沿着线𝑦=𝑥分布
   else 数据点偏离线𝑦=𝑥。

#生成1000个服从 正态分布 的数据
data_norm=stats.norm.rvs(loc=10,scale=10,size=1000)
# rvs(loc,scale,size):生成服从指定分布的随机数，
#loc：期望；scale：标准差；size：数据个数

#生成1000个服从 卡方分布 的数据
data_chi=stats.chi2.rvs(2,3,size=1000)

#画出两个概率图
fig=plt.figure(figsize=(12,6))
#正态数据
ax1=fig.add_subplot(1,2,1)
plot1=stats.probplot(data_norm,plot=ax1)
#卡方分布数据
ax2=sig.add_subplot(1,2,2)
plot2=stats.probplot(data_chi,plot=ax2)
#stats.probplot生成概率图

2.1.2两种正态性检验

正态性检验的两个假设如下：
$$H_0:\text{样本所在总体服从正态分布}\leftrightarrow H_1:\text{样本所在总体不服从正态分布}$$
⭐⭐⭐

- 小样本下使用Shapiro-Wilk检验；

无原理介绍

该检验在每一个样本值都是唯一时的检验效果最好，而一旦样本过多，难免会存在几个样本值重复的情况，这会使该方法的效用大打折扣。

样本量适用范围：样本量不得小于8，小于50最佳，小于2000效果较好，超过5000后不再适用。
⭐⭐⭐

- 大样本下使用D’Agostino’s K-squared检验

通过计算偏度（Skewness）和峰度（Kurtosis）来量化数据分布曲线与标准正态分布曲线之间的差异与不对称性

分布曲线的偏度和峰度容易受到数据量的影响，数据量越大，偏度与峰度的计算就越准确。

样本量适用范围：样本量不得小于4，除此以外越大越好。

2.1.3 同时使用多种方法判断正态性

由于数据的复杂性通常同时使用多种方法进行判断
不同方法得出的结论不同，此时就需要仔细观察数据的特征，寻找结果不一致的原因

接下来，我们在python中定义一个函数，将概率图、Shapiro-Wilk test、D’Agostino’s K-squared test结合在一起。

data_small = stats.norm.rvs(0, 1, size=30) # 小样本正态性数据集
data_large = stats.norm.rvs(0, 1, size=6000) # 大样本正态性数据集

# 定义一个正态性检验函数，它可以输出：
## 正态概率图
## 小样本Shapiro-Wilk检验的p值
## 大样本D'Agostino's K-squared检验的p值
#引入相关库
from statsmodels.stats.diagnostic import liliefors
from typing import List
#list的泛型版本
    #statsmodels.stats.diagnostic.lilliefors(x)
    #检验样本数据是否来自正态总体
  """
    输入参数
    ----------
    data : numpy数组或者pandas.Series
    show_flag : 是否显示概率图
    Returns
    -------
    两种检验的p值；概率图
    """
def check_normality(data:np.ndarray,show_flag:bool=True):
     if show_flag:
        _ = stats.probplot(data, plot=plt)
        plt.show()
    pVals = pd.Series(dtype='float64')
    # D'Agostino's K-squared test
    #大样本检验
    _, pVals['Omnibus'] = stats.normaltest(data) 

    # Shapiro-Wilk test
    #小样本检验
    _, pVals['Shapiro-Wilk'] = stats.shapiro(data)

    print(f'数据量为{len(data)}的数据集正态性假设检验的结果 : ----------------')
    print(pVals)
   
   #use 
check_normality(data_small,show_flag=True)
check_normality(data_large,show_flag=False)
 # 当样本量大于5000，会出现警告

2.2均值检验

2.2.1 单组样本均值假定的检验

检验一个样本所在总体的均值是否和一个参考值相等，这就是单组样本均值假定的检验。

Example.2 必胜中学里，陈老师班结束了一次英语考试。由于班级人数较多，短时间内很难完成批改与统计，王老师又很想知道此次班级平均分与级长定的班级均分137的目标是否有显著区别，于是他随机抽取了已经改好的10名同学的英语成绩：

136,136,134,136,131,133,142,145,137,140

问：王老师可以认为此次班级平均分与级长定的班级均分137的目标没有显著区别吗？

比较的是这个样本（10个同学的英语成绩）所代表的总体均值（班级英语成绩均值）是否与参考值137相等

two ways to solve the problem :单样本t检验与wilcoxon检验

2.2.1.1单样本t检验(正态分布)

t检验要求总体服从正态分布，即有
$$x\sim N\left(\mu ,{\sigma }^2\right)$$
要求王老师班级中每一位同学的英语成绩都服从正态分布

与z检验不同,在t检验中，总体的标准差$\sigma$是不需要事先知道的

使用p值进行假设检验的基本流程:

单样本t检验的两个假设为：
$$H_0:\mu ={\mu }_0\leftrightarrow H_1:\mu \ne {\mu }_0$$
对应的检验统计量为：
$$Teststatistics=\frac{\bar{x}-{\mu }_0}{s}$$
检验统计量服从的分布为：
$$Teststatistics\sim t_{n-1}$$
其中，$n$为样本量

2.2.1.2 wilcoxon符号秩和检验（非正态分布）

若样本数据非正态，我们应当使用wilcoxon符号秩和检验

---

秩的概念
设$x_1,\cdots ,x_n$为来自连续分布的简单随机样本，我们将它们进行从小到大排序，得到有序样本$x_{\left(1\right)}⩽\cdots ⩽x_{\left(n\right)}$。观测值$x_i$在有序样本中的序号$r$被称为$x_i$的秩。

秩和检验
设$x_1,\cdots ,x_n$为样本，对它们做绝对值变换，然后记$R_i$为$|x_i|$在$\left(|x_1|,\cdots ,|x_n|\right)$的秩。
记
$$I\left(x_i>0\right)=\{\begin{array}{ll}1,& x_i>0\\ 0,& x_i\ne 0\end{array}$$
则称
$$W^{+}=\sum _{i=1}^n{R}_{i}I\left(x_i>0\right)$$
为秩和统计量。

---

单样本均值比较的wilcoxon符号秩和检验的两个假设依旧为
$$H_0:\mu ={\mu }_0\leftrightarrow H_1:\mu \ne {\mu }_0$$

1. 计算秩和统计量
   对于待分析样本$x_1,\cdots ,x_n$，让所有样本减去对比值${\mu }_0$，得：$x_1-{\mu }_0,\cdots ,x_n-{\mu }_0$，计算出它们的秩和统计量$W^{+}$。

2. 计算检验统计量
   然后，检验统计量便可以计算为
   $$Teststatistics=\frac{W^{+}-\frac{n(n+1)}{4}}{\sqrt{\frac{n(n+1)(2n+1)}{24}}}$$

3. 检验统计量近似服从分布
   $$Teststatistics\to N\left(0,1\right)$$
   其中，$n$为样本量。

4. p值的计算与备择假设$H_1$的符号有关，这一点与t检验相同。

该方法最好在样本量大于25的情况下使用，因为这样检验统计量才近似服从正态分布

单样本均值检验代码

data=np.array([136,136,134,136,131,133,142,145,137,140])
# 定义一个单组样本均值检验函数，使它可以同时输出t检验与wilcoxon符号秩和检验的p值
 '''
    输入参数
    ----------
    data : numpy数组或者pandas.Series
    checkvalue : 想要比较的均值
    confidence : 显著性水平
    alternative : 检验类型，这取决于我们备择假设的符号:two-sided为双侧检验、greater为右侧检验、less为左侧检验

    输出
    -------
    在两种检验下的p值
    在显著性水平下是否拒绝原假设
    '''
def check_mean(data,checkvalue,confidence=0.05,alternative='two-sided'):
   pVal=pd.Series(dtype='float64')
    # 正态性数据检验-t检验
    #  ttest_1samp  测试随机样本的均值是否等于真实均值
    _,pVal['t-test']=stats.ttest_1samp(data,checkvalue,alternative=alternative)
      print('t-test------------------------')
    if pVal['t-test'] < confidence:
      print(('目标值{0:4.2f}在显著性水平{1:}下不等于样本均值(p={2:5.3f}).'.format(checkvalue,confidence,pVal['t-test'])))
    else:
      print(('目标值{0:4.2f}在显著性水平{1:}下无法拒绝等于样本均值的假设.(p={2:5.3f})'.format(checkvalue,confidence,pVal['t-test'])))

    # 非正态性数据检验-wilcoxon检验
    #stats.wilcoxon
    _, pVal['wilcoxon'] = stats.wilcoxon(data-checkvalue,alternative=alternative)
    print('wilcoxon------------------------')    
    if pVal['wilcoxon'] < confidence:
      print(('目标值{0:4.2f}在显著性水平{1:}下不等于样本均值(p={2:5.3f}).'.format(checkvalue,confidence,pVal['wilcoxon'])))
    else:
      print(('目标值{0:4.2f}在显著性水平{1:}下无法拒绝等于样本均值的假设.(p={2:5.3f})'.format(checkvalue,confidence,pVal['wilcoxon'])))
    
    return pVal

#use
check_mean(data,137,0.05)

2.2.2 两组样本的均值相等性检验

组别间独立

Example.3 陈老师在年级有一个竞争对手：王老师。王老师在他们班级在同一时间也举行了一次英语考试，且两个班用的是同一份卷子，因此两个班的英语成绩就具有比较的意义了。和陈老师一样，王老师也来不及批改和统计他们班级的英语成绩，不过他手头上也有12份已经改好的试卷，成绩分别为：

134,136,135,145,147,140,142,137,139,140,141,135

问：我们可以认为两个班级的均分是相等的吗？
一个双独立样本的均值检验
这里的独立指的是抽样意义上的独立，而不是统计意义的独立
两者在数值上是否独立（通过独立性检验判断的独立性）我们并不关心。
抽样意义上的独立：两个样本中，一个样本中的受试不能影响另一个样本中的受试。

在本例中，两个班级的授课老师不同，因此两个班学生的英语学习不会受到同一个老师的影响；两个班级考试同时进行，意味着不存在时间差让先考完的同学给后考完的同学通风报信（进而影响受试）。因此，我们可以认为这是两个独立的样本。

事实上，两个样本是否在抽样意义上独立，是没有固定答案的，因为在很多情况下，我们既不能保证两个样本间完全独立，也很难判断出两者是否存在相关性。我们只能在抽样的时候使用更加科学的抽样方法，尽可能地避免样本间的相互影响。

双样本t检验(参数检验)

两个样本的总体都服从正态分布可以使用双样本t检验

双样本t检验两个分类

• 总体标准差相等以及两种不相等
• 这两种情况下的检验统计量不相同，服从的t分布的自由度也不相同。

两种检验的假设问题是相同的：
$$H_0:{\mu }_x={\mu }_y\leftrightarrow H_1:{\mu }_x\ne {\mu }_y$$

A :若两个样本的总体方差相等

检验统计量为
$$Teststatistics=\frac{\bar{x}-\bar{y}}{s_w\sqrt{\frac{1}{m}+\frac{1}{n}}},s_w=\sqrt{\frac{1}{m+n-2}\left[\sum _{i=1}^m{\left(x_i-\bar{x}\right)}^2+\sum _{i=1}^n{\left(y_i-\bar{y}\right)}^2\right]}$$
检验统计量服从分布
$$Teststatistics\sim t_{n+m-2}$$
其中，$m$为样本$x$的样本量，$n$为样本$y$的样本量。

B :若两个样本的总体方差不相等
检验统计量为
$$Teststatistics=\frac{\bar{x}-\bar{y}}{\sqrt{\frac{s_x^2}{m}+\frac{s_y^2}{n}}}$$
检验统计量近似服从分布
$$Teststatistics\to t_l,l=\frac{{\left(\frac{s_x^2}{m}+\frac{s_y^2}{n}\right)}^2}{\frac{s_x^2}{m^2\left(m-1\right)}+\frac{s_y^2}{n^2\left(n-1\right)}}$$

判断两个总体的标准差是否相等使用方差齐性检验levene检验进行判断。

Mannwhitneyu秩和检验（非参数检验）

首先将两类样本混合在一起，对所有样本按照所考察的特征从小到大排序。
在两类样本中分别计算所得排序序号之和$T_1$和$T_2$，我们将之称作秩和。
Mannwhitneyu秩和检验的基本思想是，如果一类样本的秩和显著地比另一类小（或大），则两类样本在所考察的特征上有显著差异。

秩和检验与符号秩和检验是不一样的：前者不考虑符号的问题，而后者考虑符号问题。

def unpaired_data(group1:np.ndarray,group2:np.ndarray,confidence=0.05,alternative='two-sided'):
 """ 
    输入参数
    ----------
    group1/2 : 用于比较的两组数据
    confidence : 显著性水平
    alternative : 检验类型，这取决于我们备择假设的符号:two-sided为双侧检验、greater为右侧检验、less为左侧检验
    输出
    -------
    在两种检验下的p值
    在显著性水平下是否拒绝原假设
    """
      pVal=pd.Series(dtype='float64')

    # 先进行两组数据的**方差齐性检验**
    _,pVal['levene']=stats.levene(group1,group2)


    # t检验-若数据服从正态分布
    # 判断显著性水平与levene值大小
    if pVal['levene']<confidence:
        print('在显著性水平{0:}下，两组样本的方差不相等(p={1:.4f})，因此需要使用方差不等的t检验'.format(confidence,pVal['levene']))
        print('------------------------------------')
        _, pVal['t-test'] = stats.ttest_ind(group1, group2,equal_var=False,alternative=alternative) # 因为方差不相等，因此是False
        print('t检验p值：{}'.format(pVal['t']))
    else:
        print('在显著性水平{0:}下，不能拒绝两组样本方差相等的假设(p={1:.4f})，因此需要使用方差相等的t检验'.format(confidence,pVal['levene']))
        print('------------------------------------')
        _, pVal['t-test'] = stats.ttest_ind(group1, group2,equal_var=True,alternative=alternative) # 因为方差相等，因此是True     
        print('t检验p值：{:.3f}'.format(pVal['t-test']))   
    

    # mannwhitneyu检验-数据不服从正态检验
    _, pVal['mannwhitneyu'] = stats.mannwhitneyu(group1, group2,alternative=alternative)
    print('Mann-Whitney检验p值：{:.3f}'.format(pVal['mannwhitneyu']))
    # --- >>> STOP stats <<< ---
    
    # 两组样本均值的散点图可视化
    print('------------------------------------')
    print('两组样本均值的散点图可视化')
    plt.plot(group1, 'bx', label='group1')
    plt.plot(group2, 'ro', label='group2')
    plt.legend(loc=0)
    plt.show()

#use
# 陈老师班
group1=data
# 王老师班
group2=np.array([134,136,135,145,147,140,142,137,139,140,141,135])

unpaired_data(group1,group2)

成对组别

在进行两组间均值比较的时候，有一种特殊情况——两个样本“故意”不独立。这种情况多出现两个样本分别为同一个受试个体不同时间的受试结果
成对检验
这两个样本是“成对”的，是彼此紧密相连的。对这样两个样本进行均值比较检验，就是成对检验。

Example.4 得知了王老师班的考试均值与自己班的没有显著差异，陈老师非常生气：“我堂堂优秀教师，居然不能和隔壁老王拉开差距，岂有此理！”，于是陈老师开始了为期一周的魔鬼训练，并在一个星期后又进行了一次全班的测验。陈老师依旧抽取了上次十位同学的成绩，依次为：

139,141,137,136,135,132,141,148,145,139

问：这次班级均分与上次是否存在显著差异呢？

显然，两个样本分别为相同的同学前后两次的考试成绩，是非常典型的成对数据，因此我们可以使用成对检验。成对检验也分为两种：若总体服从正态分布，则使用成对t检验；若总体不服从正态分布，则使用成对wilcoxon秩和检验。

· 成对t检验

成对检验与独立双样本检验的假设问题是一样的
$$H_0:{\mu }_x={\mu }_y\leftrightarrow H_1:{\mu }_x\ne {\mu }_y$$

成对t检验本质上是检验${\mu }_x-{\mu }_y$的差值是否为0的单样本均值t检验

转化问题为
检验统计量就可以按照单样本均值t检验的检验统计量构造，只不过变量换成了$x-y$，目标值为0

检验统计量为：
$$Teststatistics=\frac{\bar{x}-\bar{y}}{s_d},d=x-y$$
检验统计量服从分布：
$$Teststatistics\sim t_{n-1}$$
其中，$n$为样本量。

成对wilcoxon符号秩和检验

$$H_0:{\mu }_x={\mu }_y\leftrightarrow H_1:{\mu }_x\ne {\mu }_y$$

1. 计算秩和统计量
   对于待分析样本$x_1,\cdots ,x_n$，让所有样本减去对比值${\mu }_y$，得：$x_1-{\mu }_y,\cdots ,x_n-{\mu }_y$，计算出它们的秩和统计量$W^{+}$。

2. 计算检验统计量
   然后，检验统计量便可以计算为
   $$Teststatistics=\frac{W^{+}-\frac{n(n+1)}{4}}{\sqrt{\frac{n(n+1)(2n+1)}{24}}}$$

3. 检验统计量近似服从分布
   $$Teststatistics\to N\left(0,1\right)$$
   其中，$n$为样本量。

4. p值的计算与备择假设$H_1$的符号有关，这一点与t检验相同。

该方法最好在样本量大于25的情况下使用，因为这样检验统计量才近似服从正态分布

def paired_data(group1:np.ndarray,group2:np.ndarray,confidence,alternative='two-sided'):
    """
    输入参数
    ----------
    group1/2 : 用于比较的两组数据，注意，两组数据的样本顺序必须相同
    confidence : 显著性水平
    alternative : 检验类型，这取决于我们备择假设的符号:two-sided为双侧检验、greater为右侧检验、less为左侧检验

    输出
    -------
    在两种检验下的p值
    在显著性水平下是否拒绝原假设
    """
    pVal=pd.Series(dtype='float64')

    # 配对t检验-样本服从正态分布
    _, pVal['t-test'] = stats.ttest_1samp(post - pre, 0,alternative=alternative)
    print('t-test------------------------')
    if pVal['t-test'] < confidence:
      print(('在显著性水平{0:}下，两组配对样本的均值不相等(p={1:5.3f}).'.format(confidence,pVal['t-test'])))
    else:
      print(('在显著性水平{0:}下无法拒绝等于样本均值的假设.(p={1:5.3f})'.format(confidence,pVal['t-test'])))    

    # wilcoxon秩和检验
    _, pVal['wilcoxon'] = stats.wilcoxon(group1,group2, mode='approx',alternative=alternative)
    print('wilcoxon------------------------')
    if pVal['wilcoxon'] < confidence:
      print(('在显著性水平{0:}下，两组配对样本的均值不相等(p={1:5.3f}).'.format(confidence,pVal['wilcoxon'])))
    else:
      print(('在显著性水平{0:}下无法拒绝等于样本均值的假设.(p={1:5.3f})'.format(confidence,pVal['wilcoxon'])))    
    
    return pVal
#use
# 第一次测验
pre=data
# 第二次测验
post=np.array([139,141,137,136,135,132,141,148,145,139])

paired_data(pre,post,0.05)

2.2.3 (ANOVA)方差分析-多组样本间的均值相等性检验

ANOVA(方差分析):同时对多组样本的总体进行均值检验

单因素方差分析(同时比较)

---

方差分析
单/多因素方差分析中的因素是什么意思？
举例：

A 1. 四种不同商标的同一类型产品的质量指标是否一致。A1 A2 A3 A4
B 2. 对同样的商品，采取三种不同的销售方式是否会导致显著不同的销售量。B1 B2 B3
C 3. 五个不同居民区的居民购买力是否存在显著差异。C1 C2 C3 C4 C5

---

• · 认识方差分析
  因素
  在这些例子中，商标、销售方式、居民区这些区分不同组别的依据就是因素（也叫因子），通常用大写字母$A,B,C$等表示这些因素

水平
一个因素所处的不同的状态被称为水平，用$A_1,A_2$

在这三个例子中，因素都只有一个，因此它们都是单因素方差分析；若因素有多个，则被称为多因素方差分析。

单因素方差分析而言，进行比较的样本个数，本质上就是该因素的水平个数。

方差分析的多个总体均值比较是同时比较
写成假设检验的形式有：
$$H_0:{\mu }_1={\mu }_2={\mu }_i=\cdots {\mu }_k\leftrightarrow H_1:{\mu }_i\text{不全相等}$$
其中，$k$为样本个数（因子水平个数）。

---

- 方差分析的假设

方差分析需要满足的假设：

• 每个样本的总体均为正态总体
• 各总体的标准差相同
• 从每个总体中抽取的样本相互独立

若样本满足以上条件，方差分析的结果是准确的。但是在大样本下，就算不满足假设1.与假设2.，方差分析的结果也是稳健的。

---

· 方差分析的思想

• 样本间均值的“差异程度”
• 各样本的样本内差异程度

在相同的样本间差异程度下，样本内差异程度越大，各总体间均值存在差异的可能性就越小
样本内差异程度越大，“偶然性”越大，越难以判断两个不相等的均值是否真的不相等。

举个例子：小红的考试均分是91，小刚的考试均分是89。我们假设一个非常极端的情况：他们的标准差都是0，即小红每次考试都是91，小刚每次考试都是89，那么我们似乎可以很容易地判断出，两个人的均分确实存在明显的差异；但是，如果他们的标准差都很大，高达6（方差就是36），即他们的成绩都很不稳定，这次小红考79、小刚考93，下次小红考94、小刚考70。在高达36的方差下，2分的均值差似乎没有什么说服力了。

因此，我们需要综合这两个评判指标。最简单的方法就是两者相除，即样本间均值的“差异程度”除以样本内差异程度。这就是方差分析最根本的思想。

---

· 原理解析

记$y_{ij}$为第$i$个样本的第$j$个样品，假设有r个样本，每个样本的样品数量都同为$m$，于是我们有$r\times m$个样品：
$$y_{ij},i=1,2,\cdots ,r,j=1,2,\cdots ,m,$$

既然我们要研究“差异程度”，那么观察每个样本值$y_{ij}$与平均值的偏差一定是必不可少的。数据$y_{ij}$与所有数据的总平均$\bar{y}$的偏差可用$y_{ij}-\bar{y}$表示，这个总偏差可以解构为两个偏差之和：
$$y_{ij}-\bar{y}=\left(y_{ij}-{\bar{y}}_{i\cdot }\right)+\left({\bar{y}}_{i\cdot }-\bar{y}\right)$$
其中，${\bar{y}}_{i\cdot }$为因素$A$在第$i$个水平下的样本均值（或者就是第$i$个样本的均值）。

$y_{ij}-{\bar{y}}_{i\cdot }$被称为组内偏差，它反映了第$i$个样本的组内数据与组内平均值的随机误差；还反映了第$i$个的总体均值与全总体均值的差值，它可以理解为因素$A$的第$i$个水平带来的主效应。

各$y_{ij}$间总的差异大小可以用总偏差平方和$S_T$表示
$$S_T=\sum _{i=1}^r\sum _{j=1}^m{\left(y_{ij}-\bar{y}\right)}^2,f_T=n-1,$$

仅由随机误差引起的数据间差异可用组内偏差平方和$S_e$表示
$$S_e=\sum _{i=1}^r\sum _{j=1}^m{\left(y_{ij}-{\bar{y}}_i.\right)}^2,f_e=r(m-1)=n-r$$

由于效应不同引起的数据差异可用组间偏差平方和$S_A$表示
$$S_A=m\sum _{i=1}{\left({\bar{y}}_i.-\bar{y}\right)}^2,f_A=r-1$$
其中，$f$为它们各自的自由度。

以上三种偏差平方和的大小和数据的个数（或者自由度）有关，通常而言，数据越多，偏差平方和越大，这不利于进行偏差平方和之间的比较。因此，我们需要将它们的“量纲”统一起来，方法也非常简单：让它们除以自己的自由度即可。我们定义：
$$M{S}_A=\frac{S_A}{f_A},M{S}_e=\frac{S_e}{f_e}$$
然后，我们就可以将两者相除，以判断多组样本的均值是否相等了。

单因素方差分析检验统计量为
$$Teststatistics=\frac{M{S}_A}{M{S}_e}$$
检验统计量服从分布
$$Teststatistics\sim F\left(f_A,f_e\right)$$
p值计算：(单边检验)
$$pvalue=P\left(F>Teststatistics\right)$$
尽管方差分析做的是多均值的“相等性”检验，它的p值并不是双边检验，而是单边检验。前面所说的看备择假设的符号判断p值的形式只适用于单样本与双样本均值检验。

· kruskalwallis检验——非正态性检验*(小样本,非参数)

若满足正态性假设，ANOVA的灵敏度优于kruskalwallis检验，因此我们在进行方差分析前最好做一次正态性检验

实操：对altman_910.txt数据集进行方差分析。该数据记录了3组心脏搭桥病人给予不同水平的一氧化氮通气下，他们的红细胞内叶酸水平。注：三组样本都是正态性样本。

data = np.genfromtxt('./data/altman_910.txt', delimiter=',') 
group1 = data[data[:,1]==1,0]
group2 = data[data[:,1]==2,0]
group3 = data[data[:,1]==3,0]
group1

from typing import Tuple

def anova_oneway() -> Tuple[float, float]:

    pVal=pd.Series(dtype='float64')

    # 先做方差齐性检验
    _,pVal['levene'] = stats.levene(group1, group2, group3)
    if pVal['levene']<0.05: #这里假设显著性水平为0.05
        print('警告: 方差齐性检验的p值小于0.05: p={}，方差分析结果在小样本下可能不准确'.format(pVal['levene']))
        print('-------------------------------')

    # 单因素方差分析-假设样本服从正态分布
    _, pVal['anova_oneway_normal'] = stats.f_oneway(group1, group2, group3) # 在这里输入待分析的数据
    
    print('若样本服从正态分布，单因素方差分析的p值为{}'.format(pVal['anova_oneway_normal']))
    if pVal['anova_oneway_normal'] < 0.05:
        print('检验在0.05的显著性水平下显著，多组样本中至少存在一组样本均值与其它样本的均值不相等。')

    print('---------------------------------')

    # 单因素方差分析-假设样本不服从正态分布
    _, pVal['anova_oneway_notnormal'] = stats.mstats.kruskalwallis(group1, group2, group3) # 在这里输入待分析的数据
    
    print('若样本不服从正态分布，单因素方差分析的p值为{}'.format(pVal['anova_oneway_notnormal']))
    if pVal['anova_oneway_notnormal'] < 0.05:
        print('检验在0.05的显著性水平下显著，多组样本中至少存在一组样本均值与其它样本的均值不相等。')    

    return pVal

双因素方差分析

双因素方差分析不仅仅探究两个因素是否显著，还要探究两个因素的交互项是否显著

对于这种多因素的方差分析，借用线性回归模型来解决问题可以事半功倍

· 实例讲解
我们来看一下这个测量婴儿头围hs的例子。

现在有3个婴儿的头围需要测量，我们想知道他们的头围有无显著的区别。注意，如果我们只探究这个问题，则此时就是单因素方差分析，探究的是（因素：婴儿fetus）对头围的是否有显著性影响。
然而，人工测量总是有误差的。为了保证客观性，我们请了四位观察者observer，让他们对每个婴儿的头围重复测量重复测量三次，并记录每一次测量的结果，此时，我们总共会得到334=36个样品。我们将（因素：observer）也纳入考量，同时考量两个因素以及他们的交互项的显著性。如果观察者因素没有显著差别，就说明观察者观察的结果很一致；反之，他们观察的差距非常大，可能要另外寻找其他的观察者

# statsmodel包分析的对象更多的是dataframe，前面的分析对象多为数组，这一点需要注意。
data = np.genfromtxt('./data/altman_12_6.txt', delimiter=',')
df = pd.DataFrame(data, columns=['hs', 'fetus', 'observer'])
df
import statsmodels.api as sm
from statsmodels.stats.anova import anova_lm

lm = sm.formula.ols('hs~C(fetus)+C(observer)+C(fetus):C(observer)',data=df).fit()
# 将待分析的因变量放在~左侧，因素放在~右侧
# 因素用C()括起来，交互项使用:将两者联系在一起。
anovaResults = anova_lm(lm) # 方差分析
print(anovaResults)# 3. 多元数值向量的假设检验

作业

三台车床生产同一种滚珠，我们从中各抽取13、14、16个产品，测得直径分别为：

甲车床：15.0,14.5,15.2,15.5,14.8,15.1,15.2,14.8,13.9,16.0,15.1,14.5,15.2;

乙车床：15.2,15.0,14.8,15.2,15.0,15.0,14.8,15.1,14.8,15.0,13.7,14.1,15.5,15.9;

丙车床：14.6,15.0,14.7,13.9,13.8,15.5,15.5,16.2,16.1,15.3,15.4,15,9,15.2,16.0,14.8,14.9

假设显著性水平为𝛼=0.01，问：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from typing import List
confidence=0.01
jia=np.array([15.0,14.5,15.2,15.5,14.8,15.1,15.2,14.8,13.9,16.0,15.1,14.5,15.2])
yi=np.array([15.2,15.0,14.8,15.2,15.0,15.0,14.8,15.1,14.8,15.0,13.7,14.1,15.5,15.9])
bing=np.array([14.6,15.0,14.7,13.9,13.8,15.5,15.5,16.2,16.1,15.3,15.4,15,9,15.2,16.0,14.8,14.9])

甲/乙车床生产的滚珠直径是否服从正态分布？

def check_normality(data:np.ndarray,show_flag:bool=True):
     if show_flag:
        _ = stats.probplot(data, plot=plt)
        plt.show()
    pVals = pd.Series(dtype='float64')
    # D'Agostino's K-squared test
    #大样本检验
    _, pVals['Omnibus'] = stats.normaltest(data) 

    # Shapiro-Wilk test
    #小样本检验
    _, pVals['Shapiro-Wilk'] = stats.shapiro(data)

    print(f'数据量为{len(data)}的数据集正态性假设检验的结果 : ----------------')
    print(pVals)

#use
check_normality(jia,confidence=confidence,show_flag=True)
check_normality(yi,confidence=confidence,show_flag=True)
check_normality(bing,confidence=confidence,show_flag=True)

甲/乙车床生产的滚珠直径的方差是否存在显著性差异？

    pVal=pd.Series(dtype='float64')

    # 先做方差齐性检验
    _,pVal['levene'] = stats.levene(jia,yi)
    if pVal['levene']<0.01: 
    ..
    else:
    ..

甲/乙车床生产的滚珠直径是否存在显著性差异？

def unpaired_data(group1:np.ndarray,group2:np.ndarray,confidence=0.05,alternative='two-sided'):
 """ 
    输入参数
    ----------
    group1/2 : 用于比较的两组数据
    confidence : 显著性水平
    alternative : 检验类型，这取决于我们备择假设的符号:two-sided为双侧检验、greater为右侧检验、less为左侧检验
    输出
    -------
    在两种检验下的p值
    在显著性水平下是否拒绝原假设
    """
      pVal=pd.Series(dtype='float64')

    # 先进行两组数据的**方差齐性检验**
    _,pVal['levene']=stats.levene(group1,group2)


    # t检验-若数据服从正态分布
    # 判断显著性水平与levene值大小
    if pVal['levene']<confidence:
        print('在显著性水平{0:}下，两组样本的方差不相等(p={1:.4f})，因此需要使用方差不等的t检验'.format(confidence,pVal['levene']))
        print('------------------------------------')
        _, pVal['t-test'] = stats.ttest_ind(group1, group2,equal_var=False,alternative=alternative) # 因为方差不相等，因此是False
        print('t检验p值：{}'.format(pVal['t']))
    else:
        print('在显著性水平{0:}下，不能拒绝两组样本方差相等的假设(p={1:.4f})，因此需要使用方差相等的t检验'.format(confidence,pVal['levene']))
        print('------------------------------------')
        _, pVal['t-test'] = stats.ttest_ind(group1, group2,equal_var=True,alternative=alternative) # 因为方差相等，因此是True     
        print('t检验p值：{:.3f}'.format(pVal['t-test']))   
    

    # mannwhitneyu检验-数据不服从正态检验
    _, pVal['mannwhitneyu'] = stats.mannwhitneyu(group1, group2,alternative=alternative)
    print('Mann-Whitney检验p值：{:.3f}'.format(pVal['mannwhitneyu']))
    # --- >>> STOP stats <<< ---
    
    # 两组样本均值的散点图可视化
    print('------------------------------------')
    print('两组样本均值的散点图可视化')
    plt.plot(group1, 'bx', label='group1')
    plt.plot(group2, 'ro', label='group2')
    plt.legend(loc=0)
    plt.show()

unpaired_data(jia,yi,confidence)

三台车床生产的滚珠直径是否存在显著性差异？在单因素方差分析中，本题的因素是什么？

from typing import Tuple

def anova_oneway() -> Tuple[float, float]:

    pVal=pd.Series(dtype='float64')

    # 先做方差齐性检验
    _,pVal['levene'] = stats.levene(group1, group2, group3)
    if pVal['levene']<0.05: #这里假设显著性水平为0.05
        print('警告: 方差齐性检验的p值小于0.05: p={}，方差分析结果在小样本下可能不准确'.format(pVal['levene']))
        print('-------------------------------')

    # 单因素方差分析-假设样本服从正态分布
    _, pVal['anova_oneway_normal'] = stats.f_oneway(group1, group2, group3) # 在这里输入待分析的数据
    
    print('若样本服从正态分布，单因素方差分析的p值为{}'.format(pVal['anova_oneway_normal']))
    if pVal['anova_oneway_normal'] < 0.05:
        print('检验在0.05的显著性水平下显著，多组样本中至少存在一组样本均值与其它样本的均值不相等。')

    print('---------------------------------')

    # 单因素方差分析-假设样本不服从正态分布
    _, pVal['anova_oneway_notnormal'] = stats.mstats.kruskalwallis(group1, group2, group3) # 在这里输入待分析的数据
    
    print('若样本不服从正态分布，单因素方差分析的p值为{}'.format(pVal['anova_oneway_notnormal']))
    if pVal['anova_oneway_notnormal'] < 0.05:
        print('检验在0.05的显著性水平下显著，多组样本中至少存在一组样本均值与其它样本的均值不相等。')    

    return pVal
anova_oneway(jia,yi,bing,confidence=confidence)