【直播】聊聊 AutoDev，一个可以作为 MCP 服务的 AI 辅助研发插件

接统计学整理

 置信区间：两个正态总体的情况

• 两个方差已知，求均值差的置信区间

有两个正态分布的总体X~N(\(μ_1,σ_1^2\))，Y~N(\(μ_2,σ_2^2\))，来自X的样本\(X_1,X_2,...,X_n\)，样本均值，样本方差\(S_1^2\)；来自Y的样本\(Y_1,Y_2,...,Y_n\)，样本均值，样本方差\(S_2^2\)。两个分布的置信水平都设定为1-α，\(σ_1^2\)和\(σ_2^2\)已知，我们要求的是均值差\(μ_1-μ_2\)的置信区间。

由于和也是服从正态分布的，~N(\(μ_1,{σ_1^2\over n_1}\))，~N(\(μ_2,{σ_2^2\over n_2}\))，对于两个正态分布来说，它们的和或差也是服从正态分布的，有

-~N(\(μ_1-μ_2,{σ_1^2\over n_1}+{σ_2^2\over n_2}\))

我们将该式进行平移和伸缩得到符合标准正态分布的式子

\({(mean(X)-mean(Y))-(μ_1-μ_2)\over \sqrt{{σ_1^2\over n_1}+{σ_2^2\over n_2}}}\)~N(0,1)

类比于 \({mean(X)-μ\over {σ\over \sqrt{n}}}\)~N(0,1)的置信区间\((mean(X)(+/-){σ\over \sqrt{n}}z_{α\over 2})\)，可知上式的置信区间就为

\((mean(X)-mean(Y)(+/-)\sqrt{{σ_1^2\over n_1}+{σ_2^2\over n_2}}z_{α\over 2})\)

示例4：25左右人群的月收入服从正态分布 N(\(μ_1,σ_1^2\))，35左右人群的月收入服从正态分布 N(\(μ_2,σ_2^2\))，\(σ_1,σ_2\)分别为2000和8000；我们记录了30名25岁和40名35岁个体的月收入。这30名25岁个体平均收入为16000，这40名35岁个体平均收入为25000。求\(μ_1-μ_2\)置信水平为95%的置信区间。

由题意，我们可以知道

1. 样本均值=16000，=25000
2. 样本容量\(n_1=30,n_2=40\)
3. 总体方差\(σ_1^2=2000^2,σ_2^2=8000^2\)
4. α=0.05

 这里的置信区间为

\((mean(X)-mean(Y)(+/-)\sqrt{{σ_1^2\over n_1}+{σ_2^2\over n_2}}z_{α\over 2})\)

\(=(16000-25000(+/-)\sqrt{{2000^2\over 30}+{8000^2\over 40}}*1.96)\)

\(=(-10316.56,-7683.44)\)

import numpy as np
from scipy.stats import norm
import math

def mean_diff_ci_norm_est(data1, data2, alpha, sigma1, sigma2):
    '''
    方差已知，求均值差的置信区间
    :param data1: 样本1
    :param data2: 样本2
    :param alpha: 置信水平
    :param sigma1: 标准差1
    :param sigma2: 标准差2
    :return:
    '''
    n1 = len(data1)  # 样本1的样本容量
    n2 = len(data2)  # 样本2的样本容量
    sample_mean1 = np.mean(data1)  # 样本1的样本均值
    sample_mean2 = np.mean(data2)  # 样本2的样本均值
    se = math.sqrt(sigma1**2 / n1 + sigma2**2 / n2)
    z_value = abs(norm.ppf(alpha / 2))
    return sample_mean1 - sample_mean2 - se * z_value, sample_mean1 - sample_mean2 + se * z_value

if __name__ == '__main__':

    salary_25 = np.array([1484, 785, 1598, 1366, 1716, 1020, 1716, 785, 3113, 1601])
    salary_35 = np.array([902, 4508, 3809, 3923, 4276, 2065, 1601, 553, 3345, 2182])
    print(mean_diff_ci_norm_est(salary_25, salary_35, 0.05, 2000, 8000))

运行结果

(-6308.960768849636, 3912.960768849636)

• 两个方差相等且未知，求均值差的置信区间

有两个正态分布的总体X~N(\(μ_1,σ_1^2\))，Y~N(\(μ_2,σ_2^2\))，来自X的样本\(X_1,X_2,...,X_n\)，样本均值，样本方差\(S_1^2\)；来自Y的样本\(Y_1,Y_2,...,Y_n\)，样本均值，样本方差\(S_2^2\)。两个分布的置信水平都设定为1-α，\(σ_1^2=σ_2^2=σ^2\)为未知，我们要求的是均值差\(μ_1-μ_2\)的置信区间。

由于\(σ^2\)未知，所以我们不能使用\({(mean(X)-mean(Y))-(μ_1-μ_2)\over \sqrt{{σ_1^2\over n_1}+{σ_2^2\over n_2}}}\)~N(0,1)这个统计量，只能够使用样本方差来代替总体方差，则我们可以使用

\({(mean(X)-mean(Y))-(μ_1-μ_2)\over S_W\sqrt{{1\over n_1}+{1\over n_2}}}\)~t(\(n_1+n_2-2\))

这个统计量，该统计量称为Student's t统计量，其中

\(S_W=\sqrt{(n_1-1)S_1^2+(n_2-1)S_2^2\over n_1+n_2-2}\)

它利用了两个样本的方差来得到一个新的方差，再由新的方差得到新的标准差\(S_W\)，称为合并标准差(pooled standard deviation)。

类比于 \({mean(X)-μ\over {S\over \sqrt{n}}}\)~t(n-1)的置信区间\((mean(X)(+/-){S\over \sqrt{n}}t_{α\over 2}(n-1))\)，可知这个新的统计量的置信区间就为

\((mean(X)-mean(Y)(+/-)S_W\sqrt{{1\over n_1}+{1\over n_2}}t_{α\over 2}(n_1+n_2-2))\)

示例5： 25左右人群的月收入服从正态分布 N(\(μ_1,σ_1^2\))，35左右人群的月收入服从正态分布 N(\(μ_2,σ_2^2\))，\(σ_1,σ_2\)相等但未知；我们记录了30名25岁和40名35岁个体的月收入。这30名25岁个体平均收入为16000，标准差为2500；这40名35岁个体平均收入为25000，标准差为7000。求\(μ_1-μ_2\)置信水平为95%的置信区间。

由题意，我们可以知道

1. 样本均值=16000，=25000
2. 样本容量\(n_1=30,n_2=40\)
3. 样本标准差\(S_1=2500,S_2=7000\)
4. α=0.05

 首先我们可以计算合并标准差为

\(S_W=\sqrt{(n_1-1)S_1^2+(n_2-1)S_2^2\over n_1+n_2-2}=\sqrt{(30-1)*2500^2+(40-1)*7000^2\over 30+40-2}=5546.925\)

 这里的置信区间为(\(t_{α\over 2}(30+40-2)=1.995\))

\((mean(X)-mean(Y)(+/-)S_W\sqrt{{1\over n_1}+{1\over n_2}}t_{α\over 2}(n_1+n_2-2))\)

\(=(16000-25000(+/-)5546.925*\sqrt{{1\over 30}+{1\over 40}}*1.995)\)

\(=(-11672.72,-6327.28)\)

import numpy as np
from scipy.stats import t
import math

def mean_diff_ci_t_est(data1, data2, alpha):
    '''
    方差相等未知，求均值差的置信区间
    :param data1: 样本1
    :param data2: 样本2
    :param alpha: 置信水平
    :return:
    '''
    n1 = len(data1)  # 样本1的样本容量
    n2 = len(data2)  # 样本2的样本容量
    sample_mean1 = np.mean(data1)  # 样本1的样本均值
    sample_mean2 = np.mean(data2)  # 样本2的样本均值
    s1 = np.var(data1, ddof=1)  # 样本1的方差
    s2 = np.var(data2, ddof=1)  # 样本2的方差
    sw = math.sqrt(((n1 - 1) * s1 + (n2 - 1) * s2) / (n1 + n2 - 2))
    se = sw * math.sqrt(1 / n1 + 1 / n2)
    t_value = abs(t.ppf(alpha / 2, n1 + n2 - 2))
    return sample_mean1 - sample_mean2 - se * t_value, sample_mean1 - sample_mean2 + se * t_value


if __name__ == '__main__':

    salary_25 = np.array([1484, 785, 1598, 1366, 1716, 1020, 1716, 785, 3113, 1601])
    salary_35 = np.array([902, 4508, 3809, 3923, 4276, 2065, 1601, 553, 3345, 2182])
    print(mean_diff_ci_t_est(salary_25, salary_35, 0.05))

运行结果

(-2250.697540512319, -145.30245948768106)

• 两个方差不等且未知，求均值差的置信区间

有两个正态分布的总体X~N(\(μ_1,σ_1^2\))，Y~N(\(μ_2,σ_2^2\))，来自X的样本\(X_1,X_2,...,X_n\)，样本均值，样本方差\(S_1^2\)；来自Y的样本\(Y_1,Y_2,...,Y_n\)，样本均值，样本方差\(S_2^2\)。两个分布的置信水平都设定为1-α，\(σ_1^2,σ_2^2\)为未知且不等，我们要求的是均值差\(μ_1-μ_2\)的置信区间。

这里我们将使用

\({(mean(X)-mean(Y))-(μ_1-μ_2)\over \sqrt{{S_1^2\over n_1}+{S_2^2\over n_2}}}\)~t(n)

这个统计量，该统计量称为Welch's t统计量，其中

\(n={({S_1^2\over n_1}+{S_2^2\over n_2})^2\over {({S_1^2\over n_1})^2\over n_1-1}+{({S_2^2\over n_2})^2\over n_2-1}}\)

这个新的统计量的置信区间为

\((mean(X)-mean(Y)(+/-)\sqrt{{S_1^2\over n_1}+{S_2^2\over n_2}}t_{α\over 2}(n))\)

示例6：25左右人群的月收入服从正态分布 N(\(μ_1,σ_1^2\))，35左右人群的月收入服从正态分布 N(\(μ_2,σ_2^2\))，\(σ_1,σ_2\)不等且未知；我们记录了30名25岁和40名35岁个体的月收入。这30名25岁个体平均收入为16000，标准差为2500；这40名35岁个体平均收入为25000，标准差为7000。求\(μ_1-μ_2\)置信水平为95%的置信区间。

由题意，我们可以知道

1. 样本均值=16000，=25000
2. 样本容量\(n_1=30,n_2=40\)
3. 样本标准差\(S_1=2500,S_2=7000\)
4. α=0.05

首先我们可以计算自由度为

\(n={({S_1^2\over n_1}+{S_2^2\over n_2})^2\over {({S_1^2\over n_1})^2\over n_1-1}+{({S_2^2\over n_2})^2\over n_2-1}}={({2500^2\over 30}+{7000^2\over 40})^2\over {({2500^2\over 30})^2\over 30-1}+{({7000^2\over 40})^2\over 40-1}}=51.394\)

 这里的置信区间为(\(t_{α\over 2}(51.394)=2.007\))

\((mean(X)-mean(Y)(+/-)\sqrt{{S_1^2\over n_1}+{S_2^2\over n_2}}t_{α\over 2}(n))\)

\(=(16000-25000(+/-)\sqrt{{2500^2\over 30}+{7000^2\over 40}}*2.007)\)

\(=(-11402.82,-6597.18)\)

import numpy as np
from scipy.stats import t
import math

def mean_diff_ci_t_est(data1, data2, alpha):
    '''
    方差不等未知，求均值差的置信区间
    :param data1: 样本1
    :param data2: 样本2
    :param alpha: 置信水平
    :return:
    '''
    n1 = len(data1)  # 样本1的样本容量
    n2 = len(data2)  # 样本2的样本容量
    sample_mean1 = np.mean(data1)  # 样本1的样本均值
    sample_mean2 = np.mean(data2)  # 样本2的样本均值
    s1 = np.var(data1, ddof=1)  # 样本1的方差
    s2 = np.var(data2, ddof=1)  # 样本2的方差
    se = math.sqrt(s1 / n1 + s2 / n2)
    n = (s1 / n1 + s2 / n2)**2 / ((s1 / n1)**2 / (n1 - 1) + (s2 / n2)**2 / (n2 - 1))
    t_value = abs(t.ppf(alpha / 2, n))
    return sample_mean1 - sample_mean2 - se * t_value, sample_mean1 - sample_mean2 + se * t_value


if __name__ == '__main__':

    salary_25 = np.array([1484, 785, 1598, 1366, 1716, 1020, 1716, 785, 3113, 1601])
    salary_35 = np.array([902, 4508, 3809, 3923, 4276, 2065, 1601, 553, 3345, 2182])
    print(mean_diff_ci_t_est(salary_25, salary_35, 0.05))

运行结果

(-2283.2432539051224, -112.75674609487783)

• 两个均值未知，求两个方差比的置信区间

有两个正态分布的总体X~N(\(μ_1,σ_1^2\))，Y~N(\(μ_2,σ_2^2\))，来自X的样本\(X_1,X_2,...,X_n\)，样本均值，样本方差\(S_1^2\)；来自Y的样本\(Y_1,Y_2,...,Y_n\)，样本均值，样本方差\(S_2^2\)。两个分布的置信水平都设定为1-α，\(μ_1,μ_2\)未知，我们要求的是方差比\(σ_1^2\over σ_2^2\)的置信区间。

这里我们将使用

\({{S_1^2\over S_2^2}\over {σ_1^2\over σ_2^2}}\)~F(\(n_1-1,n_2-1\))

这个统计量

\(P\{F_{1-{α\over 2}}(n_1-1,n_2-1)<{{S_1^2\over S_2^2}\over {σ_1^2\over σ_2^2}}<F_{α\over 2}(n_1-1,n_2-1)\}=1-α\)

最终我们可以得到\(σ_1^2\over σ_2^2\)的置信区间为

\(({S_1^2\over S_2^2}{1\over F_{α\over 2}(n_1-1,n_2-1)},{S_1^2\over S_2^2}{1\over F_{1-{α\over 2}}(n_1-1,n_2-1)})\)

示例7：25左右人群的月收入服从正态分布 N(\(μ_1,σ_1^2\))，35左右人群的月收入服从正态分布 N(\(μ_2,σ_2^2\))，\(μ_1,μ_2\)未知；我们记录了30名25岁和40名35岁个体的月收入。这30名25岁个体收入的标准差为2500；这40名35岁个体收入的标准差为7000。求\(σ_1^2\over σ_2^2\)置信水平为95%的置信区间。

由题意，我们可以知道

1. 样本容量\(n_1=30,n_2=40\)
2. 样本标准差\(S_1=2500,S_2=7000\)
3. α=0.05

这里的置信区间为(\(F_{0.025}(29,39)=1.962，F_{0.975}(29,39)=0.492\))

\(({S_1^2\over S_2^2}{1\over F_{α\over 2}(n_1-1,n_2-1)},{S_1^2\over S_2^2}{1\over F_{1-{α\over 2}}(n_1-1,n_2-1)})\)

\(=({2500^2\over 7000^2}*{1\over 1.962},{2500^2\over 7000^2}*{1\over 0.492})\)

\(=(0.065,0.259)\)

由该结果我们可以看出，35岁所对应的人群的工资水平变化要大于25岁人群的工资水平变化。

import numpy as np
from scipy.stats import f

def var_ratio_ci_f_est(data1, data2, alpha):
    '''
    均值未知，求方差比的置信区间
    :param data1: 样本1
    :param data2: 样本2
    :param alpha: 置信水平
    :return:
    '''
    n1 = len(data1)  # 样本1的样本容量
    n2 = len(data2)  # 样本2的样本容量
    s1 = np.var(data1, ddof=1)  # 样本1的方差
    s2 = np.var(data2, ddof=1)  # 样本2的方差
    var_ratio = s1 / s2
    f_lower_value = f.ppf(alpha / 2, n1 - 1, n2 - 1)
    f_upper_value = f.ppf(1 - alpha / 2, n1 - 1, n2 - 1)
    return var_ratio / f_upper_value, var_ratio / f_lower_value


if __name__ == '__main__':

    salary_25 = np.array([1484, 785, 1598, 1366, 1716, 1020, 1716, 785, 3113, 1601])
    salary_35 = np.array([902, 4508, 3809, 3923, 4276, 2065, 1601, 553, 3345, 2182])
    print(var_ratio_ci_f_est(salary_25, salary_35, 0.05))

运行结果

(0.05314530971751432, 0.8614126063098585)

单侧置信区间

我们之前看到的都是双侧置信区间，对于未知参数θ，我们给出两个统计量=(\(X_1,X_2,...,X_n\)) 和=(\(X_1,X_2,...,X_n\))，得到θ的双侧置信区间(,)

P{(\(X_1,X_2,...,X_n\))<θ<(\(X_1,X_2,...,X_n\))}≥1-α

在某些实际问题中，我们只关心"上限"或者"下限"，比如说电池/灯泡的平均寿命；有害物质的平均含量。

对于下限，有

P{θ> (\(X_1,X_2,...,X_n\))}≥1-α  ->  置信区间为(,+∞)

对于上限，有

P{θ< (\(X_1,X_2,...,X_n\))}≥1-α  ->  置信区间为(-∞,)

这两个就是我们所说的置信水平为1-α的单侧置信区间，为单侧置信区间的下限，为单侧置信区间的上限。

• 标准正态分布

=

在上图中，由于\(-Z_α\)左侧的面积为α，故\(-Z_α\)右侧的面积就是1-α，则\(-Z_α\)就是单侧置信区间的下限=\(-Z_α\)；同理\(Z_α\)就是单侧置信区间的上限，=\(Z_α\)。

• t分布

由上图可知=\(-t_α(n-1)\)，=\(t_α(n-1)\)。

• 卡方分布

由上图可知=\(χ_{1-α}^2(n-1)\)，=\(χ_α^2(n-1)\)。

• F分布

由上图可知=\(F_{1-α}(n_1-1,n_2-1)\)，=\(F_α(n_1-1,n_2-1)\)。

假设检验(Hypothesis Testing)

我们这里的假设检验都是在频率学派推断(Frequentist Inference)的框架下而言的，频率学派是相对于贝叶斯学派而言的，这里我们先不考虑贝叶斯学派。

假设检验：

1. 在总体的分布函数完全未知或只知其形式不知其参数的情况下，为了推断总体的某些未知特性，提出关于总体的假设。
2. 根据样本对所提出的假设做出是接受还是拒绝的决策。
3. 假设检验是做出这一决策的过程。
4. 我们不知道某一次决策是对是错，但是遵守假设检验的流程，我们至少知道犯错的频率。

假设检验的要素

• 假设(hypotheses)

1. 零假设(null hypothesis;\(H_0\))：默认为真
   1. 比如说大家长期的共识为30岁群体平均收入μ=1W
   2. 收入与性别无关\(μ_女-μ_男=0\)
2. 备择假设(alternative hypothesis;\(H_A\))：要检验的研究问题
   1. 随着时间的发展，30岁群体平均收入μ≠1W(≠,>,<)
   2.  收入与性别有关\(μ_女-μ_男≠0\)(≠,>,<)

• 检验

1. 在默认零假设为真的前提下，计算得到当前观测结果或更极端结果的概率，即P(当前观测结果或更极端结果 | \(H_0\)为真)，这是一个条件概率，有关条件概率的内容可以参考概率论整理 中的条件概率
2. 如果数据无法提供令人信服的支持备择假设的证据，则接受零假设；如果数据提供了令人信服的支持备择假设的证据，则拒绝零假设。我们可以理解成
   1.  P(当前观测结果或更极端结果 | \(H_0\)为真) > α，接受\(H_0\)
   2.  P(当前观测结果或更极端结果 | \(H_0\)为真) ≤ α，拒绝\(H_0\)
   3. α为一个阈值

频率论方法 vs 贝叶斯方法

• 贝叶斯方法

贝叶斯方法对以上三个30岁人群的平均收入不持任何偏好，然后收集数据，对这三种平均收入进行概率统计，最终得到决策，以概率最大的结果为结论。

• 频率论方法

频率论的方法会对以上两个30岁人群的平均收入存在偏好，会认为\(H_0\)是真的，然后收集数据，计算在\(H_0\)为真的条件下，有多大的可能看到当前的数据，如果这个可能性很大，我们就接受这个假设，如果这个可能性很小，我们就拒绝这个假设，最终做出决策。

正态总体均值的假设检验

• 方差已知的情况：置信区间检验法

1. \(H_0\)：μ=4    歌曲的时长服从正态分布，分布的均值μ为4
2. \(H_A\)：μ≠4    歌曲的时长服从正态分布，分布的均值μ不等于4

我们计算得到了一个置信水平为95%的歌曲平均时长的置信区间(3.8,4.2)

由于μ=4落在了这个置信区间内，所以我们接受\(H_0\)，但如果μ=3.6或者4.4，它落在了这个置信区间外，则我们就拒绝\(H_0\)。

• 方差已知的情况：p值检验法

歌曲的时长服从正态分布X~N(μ,\(σ^2\))，\(σ^2\)=1；手机里有25首歌曲，这25首歌曲的平均时长是4.5分钟

1. \(H_0\)：μ=4
2. \(H_A\)：μ>4

1.  P(当前观测结果或更极端结果 | \(H_0\)为真) > α，接受\(H_0\)
2.  P(当前观测结果或更极端结果 | \(H_0\)为真) ≤ α，拒绝\(H_0\)

这里α叫显著性水平，这里我们取α=0.05

由题意，我们可知=4.5, \(σ^2\)=1,n=25，这里我们要求的条件概率为

P(>4.5 | \(H_0\):μ=4)

由于~N(μ,\(σ^2\over n\))即为~N(4,\(1\over 25\))

则该条件概率可以转化为求上图中阴影中的面积

我们知道正态分布的概率密度函数为

\(f(x)={1\over \sqrt{2π}σ}e^{-{(x-μ)^2\over 2σ^2}}={1\over \sqrt{2π}*{1\over 5}}e^{-{(x-4)^2\over 2*{1\over 25}}}={5\over \sqrt{2π}}e^{-{25(x-4)^2\over 2}}\)

则累计概率为

P(>4.5 | \(H_0\):μ=4)=\(\int_{4.5}^{+∞}f(x)dx=\int_{4.5}^{+∞}{5\over \sqrt{2π}}e^{-{25(x-4)^2\over 2}}dx\)

这里我们使用Python来计算即可

from scipy.stats import norm

if __name__ == '__main__':

    norm_dist = norm(loc=4, scale=1 / 5)  # 创建均值为4，标准差为1/5的正态分布
    print(1 - norm_dist.cdf(4.5))  # 计算大于4.5的累积概率

运行结果

0.006209665325776159

得到的结果小于α，故\(H_0\)被拒绝。

我们也可以将该正态分布转化为标准正态分布，我们知道

\(Z={mean(X)-μ\over {σ\over \sqrt{n}}}\)~N(0,1)

由\({mean(X)-μ\over {σ\over \sqrt{n}}}={4.5-4\over {1\over \sqrt{25}}}=2.5\)

由此，我们可以将P(>4.5 | \(H_0\):μ=4)转化为P(Z>2.5)，由上图可知，当α=0.05时候的面积对应的\(Z_α=1.64\)，而2.5在1.64的右边，它所对应的面积一定小于α， 故\(H_0\)被拒绝。而P(Z>2.5)=0.006与我们用Python计算出来的结果是一致的。

我们回到的分布里，=4.5和μ=4，它们的差即距离0.5，这个距离再去除以这个分布的标准差\({σ\over \sqrt{n}}={1\over 5}\)，其实就是求距离这个对称轴有多少个标准差，如果μ=4为真，那么我们得到的一个样本均值就不应该距离μ=4太远，如果太远，我们就拒绝它。具体的标准跟α有关，每当我们确定了一个α，就可以找到一个临界值，这个临界值就是最远可以接受的样本均值距离\(H_0\)中的参数距离。

上面的这种检验我们称为右尾/边检验。

现在我们更换一下条件

1. \(H_0\)：μ=4
2. \(H_A\)：μ<4
3. α=0.05
4. =3.5
5.  \(σ^2\)=1
6. n=25

现在我们要求的条件概率是

P(<3.5 | \(H_0\):μ=4)

转化为标准正态分布

\({mean(X)-μ\over {σ\over \sqrt{n}}}={3.5-4\over {1\over \sqrt{25}}}=-2.5\)

P(<3.5 | \(H_0\):μ=4)条件概率转化为P(Z<-2.5)=0.006=p-value<α=0.05

则拒绝\(H_0\)

这样的情况，我们称为左尾/边检验；上述两种情况统称为单尾/边检验(one-tailed test,one-sided test)

我们再更换一下条件

1. \(H_0\)：μ=4
2. \(H_A\)：μ≠4
3. α=0.05
4. =4.5或=3.5(这是两个样本)
5.  \(σ^2\)=1
6. n=25

现在我们要求的条件概率是

P(>4.5 or <3.5 | \(H_0\):μ=4)

转化为标准正态分布

1. \({mean(X)-μ\over {σ\over \sqrt{n}}}={4.5-4\over {1\over \sqrt{25}}}=2.5\)
2. \({mean(X)-μ\over {σ\over \sqrt{n}}}={3.5-4\over {1\over \sqrt{25}}}=-2.5\)

P(>4.5 or <3.5 | \(H_0\):μ=4)条件概率转化为

P(Z>2.5 or Z<-2.5)=0.006+0.006=0.012=p-value<α=0.05， 拒绝\(H_0\)

这种使用Z分数的检验称为Z检验，这种情况称为双尾/边检验(two-tailed test,two-sided test)。

import numpy as np
import math
from scipy.stats import norm

def z_test(data, tail, mu, sigma):
    '''
    z检验(方差已知)
    :param data: 样本
    :param mu: 总体均值
    :param sigma: 总体标准差
    :return:
    '''
    assert tail in ['both', 'left', 'right'], 'tail should be one of "both", "left", "right"'
    mean_val = np.mean(data)  # 样本均值
    n = len(data)  # 样本容量
    se = sigma / math.sqrt(n)  # 标准误
    z_value = (mean_val - mu) / se
    if tail == 'both':
        p_value = 2 * (1 - norm.cdf(abs(z_value)))
    elif tail == 'left':
        p_value = norm.cdf(z_value)
    else:
        p_value = 1 - norm.cdf(z_value)
    return z_value, p_value

if __name__ == '__main__':

    data = np.array([41, 36, 12, 18, 23, 19, 8, 16, 11, 14, 18, 14, 34, 6, 30, 11, 1, 11, 4, 32])
    print(z_test(data, tail='both', mu=35, sigma=5))

运行结果

(-15.249983606548566, 0.0)

由结果可知，p-value=0< α=0.05，拒绝\(H_0\)。

• 方差未知的情况：p 值检验法

歌曲的时长服从正态分布 X~N(μ,\(σ^2\))，\(σ^2\)未知；手机里有25首歌曲，这25首歌曲的平均时长是4.5分钟，标准差为1分钟

1.  \(H_0\)：μ=4
2. \(H_A\)：μ>4
3. α=0.05
4. =4.5
5.  S=1
6. n=25

之前我们说过这个统计量

\(T={mean(X)-μ\over {S\over \sqrt{n}}}\)~t(n-1)

我们将转换为T值

\({mean(X)-μ\over {S\over \sqrt{n}}}={4.5-4\over {1\over \sqrt{25}}}=2.5\)

P(>4.5 | \(H_0\):μ=4)条件概率转化为P(T>2.5)

α=0.05面积所对应的T值为1.71，因为2.5>1.71，它所对应的面积P(T>2.5)=0.01<α=0.05，故拒绝\(H_0\)。

这里跟之前的不同在于之前是利用了标准正态分布来进行检验，而这里是利用了t分布来进行检验，它们的流程是一样的，这样的检验称为t检验，对于这个单个总体均值的情况，该t检验又称为单样本t检验(one-sample t-test)。

现在我们更换一下条件

1.  \(H_0\)：μ=4
2. \(H_A\)：μ>4
3. α=0.05
4. =4.1
5.  S=1
6. n=25

\({mean(X)-μ\over {S\over \sqrt{n}}}={4.1-4\over {1\over \sqrt{25}}}=0.5\)

  P(>4.1 | \(H_0\):μ=4)条件概率转化为P(T>0.5)

α=0.05面积所对应的T值为1.71，因为0.5<1.71，它所对应的面积P(T>0.5)=0.31>α=0.05，故接受\(H_0\)。

import numpy as np
import math
from scipy.stats import t

def t_test(data, tail, mu):
    '''
    t检验(方差未知)
    :param data: 样本
    :param mu: 总体均值
    :return:
    '''
    assert tail in ['both', 'left', 'right'], 'tail should be one of "both", "left", "right"'
    mean_val = np.mean(data)  # 样本均值
    n = len(data)  # 样本容量
    df = n - 1  # t分布的自由度
    s = np.std(data, ddof=1)  # 样本标准差
    se = s / math.sqrt(n)  # 标准误
    t_value = (mean_val - mu) / se
    if tail == 'both':
        p_value = 2 * (1 - t.cdf(abs(t_value), df))
    elif tail == 'left':
        p_value = t.cdf(t_value, df)
    else:
        p_value = 1 - t.cdf(t_value, df)
    return t_value, p_value

if __name__ == '__main__':

    data = np.array([41, 36, 12, 18, 23, 19, 8, 16, 11, 14, 18, 14, 34, 6, 30, 11, 1, 11, 4, 32])
    print(t_test(data, tail='both', mu=35))

运行结果

(-6.752179134841297, 1.886684948049222e-06)

由结果可知，p-value=1.886684948049222e-06< α=0.05，拒绝\(H_0\)。

两个总体均值差的检验

• 方差已知的情况：p值检验法

25岁左右人群的月收入服从正态分布N(\(μ_1,σ_1^2\))，35岁左右人群的月收入服从正态分布N(\(μ_2,σ_2^2\))，\(σ_1,σ_2\)分别为2000和8000；我们记录了30名25岁和40名35岁个体的月收入。这30名25岁个体平均月收入为16000，这40名35岁个体平均收入为25000。

1. \(H_0:μ_1-μ_2=0\)
2. \(H_A:μ_1-μ_2≠0\)
3. α=0.05
4. =16000,=25000
5. \(σ_1=2000,σ_2=8000\)
6. \(n_1=30,n_2=40\)

P(当前观测结果或更极端结果 | \(H_0\)为真)

即

P(-<-9000 or ->9000 | \(H_0:μ_1-μ_2=0\))

这里的-9000为16000-25000得来，因为是双尾检验，所以还要检查右尾->9000。

-~N(\(μ_1-μ_2,{σ_1^2\over n_1}+{σ_2^2\over n_2}\))

转化为标准正态分布为

\({(mean(X)-mean(Y))-(μ_1-μ_2)\over \sqrt{{σ_1^2\over n_1}+{σ_2^2\over n_2}}}\)~N(0,1)

\({(mean(X)-mean(Y))-(μ_1-μ_2)\over \sqrt{{σ_1^2\over n_1}+{σ_2^2\over n_2}}}={-9000-0\over \sqrt{{2000^2\over 30}+{8000^2\over 40}}}=-6.84\)

再将9000代入为

\({(mean(X)-mean(Y))-(μ_1-μ_2)\over \sqrt{{σ_1^2\over n_1}+{σ_2^2\over n_2}}}={9000-0\over \sqrt{{2000^2\over 30}+{8000^2\over 40}}}=6.84\)

P(-<-9000 or ->9000 | \(H_0:μ_1-μ_2=0\))可以转化为P(Z>6.84 or Z<-6.84)

P(Z>6.84 or Z<-6.84)=3.96e-12 * 2=7.92e-12 < 0.05，故拒绝\(H_0\)。

import numpy as np
import math
from scipy.stats import norm

def z_test(data1, data2, tail, mu=0, sigma1=1, sigma2=1):
    '''
    z检验：方差已知
    :param data1: 样本1
    :param data2: 样本2
    :param tail: 双尾，左尾，右尾
    :param mu: 总体均值差
    :param sigma1: 样本1总体标准差
    :param sigma2: 样本2总体标准差
    :return:
    '''
    assert tail in ['both', 'left', 'right'], 'tail should be one of "both", "left", "right"'
    mean_val1 = np.mean(data1)  # 样本1均值
    mean_val2 = np.mean(data2)  # 样本2均值
    mean_diff = mean_val1 - mean_val2
    n1 = len(data1)  # 样本1容量
    n2 = len(data1)  # 样本1容量
    se = math.sqrt(sigma1**2 / n1 + sigma2**2 / n2)  # 标准误
    z_value = (mean_diff - mu) / se
    if tail == 'both':
        p_value = 2 * (1 - norm.cdf(abs(z_value)))
    elif tail == 'left':
        p_value = norm.cdf(z_value)
    else:
        p_value = 1 - norm.cdf(z_value)
    return z_value, p_value

if __name__ == '__main__':

    data1 = np.array([41, 36, 12, 18, 23, 19, 8, 16, 11, 14, 18, 14, 34, 6, 30, 11, 1, 11, 4, 32])
    data2 = np.array([23, 45, 115, 37, 29, 71, 39, 23, 21, 37, 20, 12, 13, 135, 49, 32, 64, 40, 77, 97])
    print(z_test(data1, data2, tail='both', mu=0, sigma1=5, sigma2=15))

运行结果

(-8.76812408671319, 0.0)

 由结果可知，p-value=0< α=0.05，拒绝\(H_0\)。

• 方差未知且相等的情况：p值检验法

25岁左右人群的月收入服从正态分布N(\(μ_1,σ_1^2\))，35岁左右人群的月收入服从正态分布N(\(μ_2,σ_2^2\))，\(σ_1,σ_2\)相等且未知；我们记录了30名25岁和40名35岁个体的月收入。这30名25岁个体平均月收入为16000，标准差为2500，这40名35岁个体平均收入为25000，标准差为7000.

1. \(H_0:μ_1-μ_2=0\)
2. \(H_A:μ_1-μ_2≠0\)
3. α=0.05
4. =16000,=25000
5. \(S_1=2500,S_2=7000\)
6. \(n_1=30,n_2=40\)

P(当前观测结果或更极端结果 | \(H_0\)为真)

依然为

P(-<-9000 or ->9000 | \(H_0:μ_1-μ_2=0\))

\({(mean(X)-mean(Y))-(μ_1-μ_2)\over S_W\sqrt{{1\over n_1}+{1\over n_2}}}\)~t(\(n_1+n_2-2)\)

\(S_W=\sqrt{(n_1-1)S_1^2+(n_2-1)S_2^2\over n_1+n_2-2}=\sqrt{(30-1)*2500^2+(40-1)*7000^2\over 30+40-2}=5547\)

\({(mean(X)-mean(Y))-(μ_1-μ_2)\over S_W\sqrt{{1\over n_1}+{1\over n_2}}}={(+/-)9000-0\over 5547*\sqrt{{1\over 30}+{1\over 40}}}=(+/-)6.72\)

P(-<-9000 or ->9000 | \(H_0:μ_1-μ_2=0\))可以转化为P(T>6.72 or T<-6.72)

P(T>6.72 or T<-6.72)=2.26e-9*2=4.52e-9<0.05 ，故拒绝\(H_0\)。由于是t分布，所以这里也是t检验，该t检验又称为两个独立样本的t检验(Two independent sample t-test)。

import numpy as np
import math
from scipy.stats import t

def t_test(data1, data2, tail, mu=0):
    '''
    t检验(方差未知且相等)
    :param data1: 样本1
    :param data2: 样本2
    :param tail: 双尾，左尾，右尾
    :param mu: 总体均值差
    :return:
    '''
    assert tail in ['both', 'left', 'right'], 'tail should be one of "both", "left", "right"'
    mean_val1 = np.mean(data1)  # 样本1均值
    mean_val2 = np.mean(data2)  # 样本2均值
    mean_diff = mean_val1 - mean_val2
    n1 = len(data1)  # 样本1容量
    n2 = len(data2)  # 样本2容量
    df = n1 + n2 - 2  # t分布的自由度
    s1 = np.var(data1, ddof=1)  # 样本1方差
    s2 = np.var(data2, ddof=1)  # 样本2方差
    sw = math.sqrt(((n1 - 1) * s1 + (n2 - 1) * s2) / df)
    se = sw * math.sqrt(1 / n1 + 1 / n2)  # 标准误
    t_value = (mean_diff - mu) / se
    if tail == 'both':
        p_value = 2 * (1 - t.cdf(abs(t_value), df))
    elif tail == 'left':
        p_value = t.cdf(t_value, df)
    else:
        p_value = 1 - t.cdf(t_value, df)
    return t_value, p_value

if __name__ == '__main__':

    data1 = np.array([41, 36, 12, 18, 23, 19, 8, 16, 11, 14, 18, 14, 34, 6, 30, 11, 1, 11, 4, 32])
    data2 = np.array([23, 45, 115, 37, 29, 71, 39, 23, 21, 37, 20, 12, 13, 135, 49, 32, 64, 40, 77, 97])
    print(t_test(data1, data2, tail='both', mu=0))

运行结果

(-3.8382532331766246, 0.0004548911808865963)

  由结果可知，p-value=0.0004548911808865963< α=0.05，拒绝\(H_0\)。

• 方差未知且不等的情况：p值检验法

25岁左右人群的月收入服从正态分布N(\(μ_1,σ_1^2\))，35岁左右人群的月收入服从正态分布N(\(μ_2,σ_2^2\))，\(σ_1,σ_2\)不等且未知；我们记录了30名25岁和40名35岁个体的月收入。这30名25岁个体平均月收入为16000，标准差为2500，这40名35岁个体平均收入为25000，标准差为7000.

1. \(H_0:μ_1-μ_2=0\)
2. \(H_A:μ_1-μ_2≠0\)
3. α=0.05
4. =16000,=25000
5. \(S_1=2500,S_2=7000\)
6. \(n_1=30,n_2=40\)

P(当前观测结果或更极端结果 | \(H_0\)为真)

依然为

P(-<-9000 or ->9000 | \(H_0:μ_1-μ_2=0\))

\({(mean(X)-mean(Y))-(μ_1-μ_2)\over \sqrt{{S_1^2\over n_1}+{S_2^2\over n_2}}}\)~t(n)

\(n={({S_1^2\over n_1}+{S_2^2\over n_2})^2\over {({S_1^2\over n_1})^2\over n_1-1}+{({S_2^2\over n_2})^2\over n_2-1}}={({2500^2\over 30}+{7000^2\over 40})^2\over {({2500^2\over 30})^2\over 30-1}+{({7000^2\over 40})^2\over 40-1}}=51.394\)

\({(mean(X)-mean(Y))-(μ_1-μ_2)\over \sqrt{{S_1^2\over n_1}+{S_2^2\over n_2}}}={(+/-)9000-0\over \sqrt{{2500^2\over 30}+{7000^2\over 40}}}=(+/-)7.52\)

P(-<-9000 or ->9000 | \(H_0:μ_1-μ_2=0\))可以转化为P(T>7.52 or T<-7.52)

P(T>7.52 or T<-7.52)=3.93e-10*2=7.86e-10<0.05，故拒绝\(H_0\)。

import numpy as np
import math
from scipy.stats import t

def t_test(data1, data2, tail, mu=0):
    '''
    t检验(方差未知且不等)
    :param data1: 样本1
    :param data2: 样本2
    :param tail: 双尾，左尾，右尾
    :param mu: 总体均值差
    :return:
    '''
    assert tail in ['both', 'left', 'right'], 'tail should be one of "both", "left", "right"'
    mean_val1 = np.mean(data1)  # 样本1均值
    mean_val2 = np.mean(data2)  # 样本2均值
    mean_diff = mean_val1 - mean_val2
    n1 = len(data1)  # 样本1容量
    n2 = len(data2)  # 样本2容量
    s1 = np.var(data1, ddof=1)  # 样本1方差
    s2 = np.var(data2, ddof=1)  # 样本2方差
    df = (s1 / n1 + s2 / n2) ** 2 / ((s1 / n1) ** 2 / (n1 - 1) + (s2 / n2) ** 2 / (n2 - 1))  # t分布的自由度
    se = math.sqrt(s1 / n1 + s2 / n2)  # 标准误
    t_value = (mean_diff - mu) / se
    if tail == 'both':
        p_value = 2 * (1 - t.cdf(abs(t_value), df))
    elif tail == 'left':
        p_value = t.cdf(t_value, df)
    else:
        p_value = 1 - t.cdf(t_value, df)
    return t_value, p_value

if __name__ == '__main__':

    data1 = np.array([41, 36, 12, 18, 23, 19, 8, 16, 11, 14, 18, 14, 34, 6, 30, 11, 1, 11, 4, 32])
    data2 = np.array([23, 45, 115, 37, 29, 71, 39, 23, 21, 37, 20, 12, 13, 135, 49, 32, 64, 40, 77, 97])
    print(t_test(data1, data2, tail='both', mu=0))

运行结果

(-3.8382532331766246, 0.000836440983258413)

  由结果可知，p-value=0.000836440983258413< α=0.05，拒绝\(H_0\)。

• 基于成对数据的检验

在我们之前所介绍的对两个样本进行检验中，都有一个前提条件就是这两个样本是相互独立的。所谓相互独立是指我们不能基于第一个样本获得第二个样本的任何信息，在上表中，显然训练前和训练后的数学成绩不满足相互独立这个条件，训练后的数学成绩肯定是围绕着训练前的数学成绩浮动的，比如第一个人训练前得70分，我们大概率不会猜他训练后得100分，而一个人训练前得95分，我们可以猜他训练后得100分，即便他考的不理想，我们也会猜他的成绩比较接近95，而不是接近70。这样的数据就叫做成对数据。

对于这种不满足相互独立的成对数据，我们可以使用差值来作为单独的样本，去检验该变量的总体。那么该总体的均值就为

1. 参数：所有参加训练的学生的两次成绩之差的平均值\(μ_{post-pre}\)
2. 点估计：样本包含的学生的两次成绩之差的平均值
3. \(H_0:μ_{post-pre}=0\)  训练前后数学成绩没有差别
4. \(H_A:μ_{post-pre}≠0\)  训练前后数学成绩有差别
5. \(H_A:μ_{post-pre}>0\)  训练前后数学成绩有提高
6. \(H_A:μ_{post-pre}<0\)  训练前后数学成绩有下降

从所有参加训练的学生中抽取了50人，采集了训练前后的数学成绩；对每一个学生计算数学成绩之差(训练后-训练前)后，得到成绩之差的均值为10，标准差为3；问题：训练是否有效提高了数学成绩？

1. \(H_0:μ_{post-pre}=0\)
2. \(H_A:μ_{post-pre}>0\)
3. α=0.05
4. =10
5. S=3
6. n=50

P(当前观测结果或更极端结果 | \(H_0\)为真)

即为

P(>10 | \(H_0:μ_{post-pre}=0\))

这是一个总体方差未知的情况，则我们计算T值

\({mean(X)-μ\over {S\over \sqrt{n}}}={10-0\over {3\over \sqrt{50}}}=23.57\)

P(T>23.57)=1.09e-28<0.05，拒绝\(H_0\)

这个检验也是t检验，该t检验称为配对t检验(paired t-test)。

• 单个 vs 两个总体均值的检验

单个总体均值的检验以及配对t检验     <-->    1个数值变量

上表中身高就是单个总体，身高差为配对t检验

两个总体均值的检验    <-->    1个数值变量和1个分类变量

上表中身高就是数值变量，性别就是分类变量，依据分类变量，我们可以将身高变量分成两组，一组是男生的身高数据，一组是女生的身高数据。这里我们真正感兴趣的是身高这个变量，为响应变量(Response variable)；而性别为解释变量(Explanatory variable)，我们引入性别这个因素是为了看能否解释身高的变化。

import numpy as np
import math
from scipy.stats import t

def t_test_paired(data1, data2, tail, mu=0):
    '''
    t检验(成对数据检验)
    :param data1: 样本1
    :param data2: 样本2
    :param tail: 双尾，左尾，右尾
    :param mu: 总体均值差
    :return:
    '''
    data = np.array([e1 - e2 for (e1, e2) in zip(data1, data2)])
    assert tail in ['both', 'left', 'right'], 'tail should be one of "both", "left", "right"'
    mean_val = np.mean(data)  # 样本均值
    n = len(data)  # 样本容量
    df = n - 1  # t分布的自由度
    s = np.std(data, ddof=1)  # 样本标准差
    se = s / math.sqrt(n)  # 标准误
    t_value = (mean_val - mu) / se
    if tail == 'both':
        p_value = 2 * (1 - t.cdf(abs(t_value), df))
    elif tail == 'left':
        p_value = t.cdf(t_value, df)
    else:
        p_value = 1 - t.cdf(t_value, df)
    return t_value, p_value

if __name__ == '__main__':

    data1 = np.array([41, 36, 12, 18, 23, 19, 8, 16, 11, 14, 18, 14, 34, 6, 30, 11, 1, 11, 4, 32])
    data2 = np.array([23, 45, 115, 37, 29, 71, 39, 23, 21, 37, 20, 12, 13, 135, 49, 32, 64, 40, 77, 97])
    print(t_test_paired(data1, data2, tail='both', mu=0))

运行结果

(-3.561497630199985, 0.0020830862083978108)

 由结果可知，p-value=0.0020830862083978108< α=0.05，拒绝\(H_0\)。

置信区间与假设检验的关系

显著性水平为α的双尾假设检验   <=>   置信水平为1-α的双侧置信区间

这两者的关系是等价的，比如说

1. \(H_0:μ=μ_0\)
2. \(H_A:μ≠μ_0\)
3. α=0.05

上图是一个标准正态分布的图形，由于是双尾检验，α=0.05，所以我们分摊到两边的面积就是各0.025，它们对应的Z值就是-1.96和1.96。由于样本均值

 ~N(μ,\(σ^2\over n\))

这个正态分布，经过平移和伸缩，我们将该正态分布转换为标准正态分布，就有

\(Z={mean(X)-μ\over {σ\over \sqrt{n}}}\)~N(0,1)

通过计算这个Z值，我们只需要将该值与-1.96和1.96做比较，如果这个Z值落在了-1.96左侧或者1.96右侧，我们就拒绝零假设\(H_0\)，否则我们就接受\(H_0\)。

这个是p检验的过程，而双侧置信区间为

\((mean(X)(+/-){σ\over \sqrt{n}}z_{α\over 2})\)

用图形来表示即为

如果\(μ_0\)处于区间内就接受\(H_0\)，否则就拒绝\(H_0\)。

由\(Z={mean(X)-μ_0\over {σ\over \sqrt{n}}}\)，可得

\(μ_0=mean(X)-Z{σ\over \sqrt{n}}\)

当Z≥1.96，\(μ_0≤mean(X)-1.96{σ\over \sqrt{n}}\)

从p检验的标准正态分布图来看，\(μ_0\)进入了1.96右侧的阴影区域，表示该拒绝\(H_0\)；从置信区间的图来看，\(μ_0\)落入了置信区间下限的左侧，所以同样拒绝\(H_0\)。所以两者是一致的。

 当Z≤-1.96，\(μ_0≥mean(X)+1.96{σ\over \sqrt{n}}\)

从p检验的标准正态分布图来看，\(μ_0\)进入了-1.96左侧的阴影区域，表示该拒绝\(H_0\)；从置信区间的图来看，\(μ_0\)落入了置信区间上限的右侧，所以同样拒绝\(H_0\)。所以这两者也是一致的。

我们可以看到，p检验的标准正态分布图左尾的拒绝域与置信区间右侧的拒绝区间联系了起来；p检验的标准正态分布图右尾的拒绝域与置信区间左侧的拒绝区间联系了起来；p检验的标准正态分布图中间部分的接受域与置信区间的接受区间联系了起来。

所以我们可以看到基于置信区间的假设检验与基于p检验的假设检验的本质是一样的。

单尾和单侧区间与上面是差不多的，这里就不写了。最终可得

如果假设检验的结果是拒绝\(H_0\)，则置信区间应该不包含\(H_0\)中的参数值；如果假设检验的结果是接受\(H_0\)，则置信区间应该包含\(H_0\)中的参数值。

正态总体方差的假设检验

• 单个总体方差的检验

25岁左右人群的月收入在3年前服从标准差为2(千元)的正态分布；我们现收集了30名25岁个体的月收入，标准差为2.5(千元)。问能够推断现在25岁人群月收入的波动性与3年前有显著变化吗？

1. \(H_0:σ^2=4\)
2. \(H_A:σ^2≠4\)
3. α=0.05
4. S=2.5
5. n=30

P(当前观测结果或更极端结果 | \(H_0\)为真)

即为

P(\(S^2>4\) or \(S^2<-4\) | \(H_0:σ^2=4\))

由 \({(n-1)S^2\over σ^2}\)~\(χ^2(n-1)\)计算卡方值为

\({(n-1)S^2\over σ^2}={(30-1)*2.5^2\over 4}=45.31\)

转化为卡方分布为

p-value=2*P(\(χ^2>45.31\))

由于卡方概率密度函数并不是对称的，所以左尾面积为0.025和右尾面积为0.025的卡方值并不相等，左尾为16.05，右尾为45.72。在右尾，我们计算出来的卡方值45.31是小于45.72，且相同的面积在左尾的卡方值为16.25，它是大于16.05的。

 2*P(\(χ^2>45.31\))=2*0.027=0.054>0.05，接受\(H_0\)。该检验称为卡方检验。

import numpy as np
from scipy.stats import chi2

def chi2_test(data, tail, sigma2):
    '''
    卡方检验：单个总体方差的检验
    :param data: 样本
    :param tail: 双尾，左尾，右尾
    :param sigma2: 总体方差
    :return:
    '''
    assert tail in ['both', 'left', 'right'], 'tail should be one of "both", "left", "right"'
    n = len(data)  # 样本容量
    df = n - 1  # 卡方分布的自由度
    s = np.var(data, ddof=1)  # 样本方差
    chi2_value = (n - 1) * s / sigma2
    if tail == 'both':
        p_value = 2 * min(1 - chi2.cdf(chi2_value, df), chi2.cdf(chi2_value, df))
    elif tail == 'left':
        p_value = chi2.cdf(chi2_value, df)
    else:
        p_value = 1 - chi2.cdf(chi2_value, df)
    return chi2_value, p_value

if __name__ == '__main__':

    data = np.array([41, 36, 12, 18, 23, 19, 8, 16, 11, 14, 18, 14, 34, 6, 30, 11, 1, 11, 4, 32])
    print(chi2_test(data, tail='both', sigma2=5))

运行结果

(484.59000000000003, 0.0)

 由结果可知，p-value=0.0< α=0.05，拒绝\(H_0\)。

• 两个总体方差比的检验

我们现在收集了30名25岁个体的月收入，标准差为2.5(千元)；收集了40名35岁个体的月收入，标准差为4.2(千元)。问能够推断25岁人群月收入的波动小于35岁人群月收入的波动吗？

1. \(H_0:{σ_1^2\over σ_2^2}=1\)
2. \(H_A:{σ_1^2\over σ_2^2}<1\)
3. α=0.05
4. \(S_1=2.5,S_2=4.2\)
5. \(n_1=30,n_2=40\)

P(当前观测结果或更极端结果 | \(H_0\)为真)

即为

P(\({S_1^2\over S_2^2}<1\) | \(H_0:{σ_1^2\over σ_2^2}=1\))

\({{S_1^2\over S_2^2}\over {σ_1^2\over σ_2^2}}\)~F(\(n_1-1,n_2-1\)

\({{S_1^2\over S_2^2}\over {σ_1^2\over σ_2^2}}={{2.5^2\over 4.2^2}\over 1}=0.35\)

转化为F分布为P(F<0.35)

由于这只是一个左尾检测，面积为0.05在左尾的卡方值为0.55，我们计算出来的卡方值0.35落在了0.55的左侧。p-value=P(F<0.35)=0.002<0.05，拒绝\(H_0\)。这也是一个卡方检验。

有的时候，我们可以使用该检验来检查是否满足方差相等这个前提条件，如我们想检验25岁群体和35岁群体的平均月收入是否有显著差异，但是我们不知道这两个群体的方差是否相等(方差是否相等决定了应该使用哪种t检验)，所以我们先对方差是否相等进行检验。两总体的方差相等也称为两总体具有方差齐性。

import numpy as np
from scipy.stats import f

def f_test(data1, data2, tail, ratio=1):
    '''
    F检验(两个总体方差比的检验)
    :param data1: 样本1
    :param data2: 样本2
    :param tail: 双尾，左尾，右尾
    :param ratio: 总体方差比
    :return:
    '''
    assert tail in ['both', 'left', 'right'], 'tail should be one of "both", "left", "right"'
    n1 = len(data1)  # 样本1的容量
    n2 = len(data2)  # 样本2的容量
    df1 = n1 - 1  # 第一个样本的F自由度
    df2 = n2 - 1  # 第二个样本的F自由度
    s1 = np.var(data1, ddof=1)  # 样本1的方差
    s2 = np.var(data2, ddof=1)  # 样本2的方差
    f_value = s1 / s2 / ratio  # F值
    if tail == 'both':
        p_value = 2 * min(1 - f.cdf(f_value, df1, df2), f.cdf(f_value, df1, df2))
    elif tail == 'left':
        p_value = f.cdf(f_value, df1, df2)
    else:
        p_value = 1 - f.cdf(f_value, df1, df2)
    return f_value, p_value

if __name__ == '__main__':

    data1 = np.array([41, 36, 12, 18, 23, 19, 8, 16, 11, 14, 18, 14, 34, 6, 30, 11, 1, 11, 4, 32])
    data2 = np.array([23, 45, 115, 37, 29, 71, 39, 23, 21, 37, 20, 12, 13, 135, 49, 32, 64, 40, 77, 97])
    print(f_test(data1, data2, tail='both', ratio=1))

运行结果

(0.10833692898933374, 1.1011019535150277e-05)

 由结果可知，p-value=1.1011019535150277e-05< α=0.05，拒绝\(H_0\)。

决策错误与统计功效

P(当前观测结果或更极端结果 | \(H_0\)为真)

1. p>α，接受\(H_0\)(结果不显著)
2. p≤α，拒绝\(H_0\)(结果显著)

上表中一类错误为当\(H_0\)为真时，我们却拒绝了\(H_0\)，它的概率为P(reject \(H_0\) | \(H_0\) true)=α；二类错误为当\(H_A\)为真时，我们却接受了\(H_0\)，它的概率为P(accept \(H_0\) | \(H_A\) true)=β。α和β两者是相互制约的。

我们来看一个例子

1. \(H_0\)：被告无罪
2. \(H_A\)：被告有罪

它可能产生的决策错误有

1. 被告是无辜的，但是法院判决被告有罪；     一类错误
2. 被告不是无辜的，但是法院判决被告无罪。     二类错误

由于一类错误和二类错误是相互制约的，减小其中一种错误的概率就会增加其中另一种错误的概率。我们在面对有可能犯错的情况下作出合理的决策主要取决于哪类错误的后果更严重。

1. 如果一类错误更严重，则选择低的显著性水平(例如α=0.01)。在这个示例中，一类错误明显要比二类错误更加严重，所以我们可以将α选的更小一些，让我们犯一类错误的概率更低。
2. 如果二类错误更严重，则选择高的显著性水平(例如α=0.10)。例如病人是否有癌症，一类错误是病人没有癌症，但医生判定病人有癌症；二类错误是病人有癌症，但医生判定病人没有癌症。在这两种错误中，二类错误更加严重。

这里的1-β是一个重要的概念，叫做统计功效(statistical power)。我们希望犯一类错误和二类错误的概率都尽量小，也就是我们想同时保证α和β都小。我们对犯一类错误的控制是直接的，也就是在做统计推断的时候，就已经设计好了α的取值。如果想让一类错误概率尽量小的话，只需要让α尽量小就可以了。

我们知道α和β之间有相互制约的关系，减小α在一定程度上就会增大β，我们想减小β，但是我们却不直接去讨论β，我们要讨论的是统计功效，这个统计功效越大的话，就意味着β越小，也就是说犯二类错误概率的降低=统计功效的升高。

• 统计功效

如果备择假设是真的，犯二类错误的概率是多少？

1. \(H_A\)为真时，如果真值与零假设参数值非常接近，则很难检测到差异(即更难拒绝\(H_0\)，更容易犯错误)；如果真值与零假设参数值非常远，则很容易检测到差异(即更易拒绝\(H_0\)，更不容易犯错误)。
2. Power=P(拒绝\(H_0\) ｜ \(H_A\)为真)的大小取决于真值与零假设参数值的距离

以往的资料告诉我们，一班和二班(各50人)数学成绩的标准差均为8；一班接受数学训练，二班没有接受数学训练，现收集了两班的数学成绩。问两班成绩之差至少为多大可以让我们在显著性水平为0.05时拒绝零假设？

1. \(H_0:μ_1-μ_2=0\)
2. \(H_A:μ_1-μ_2≠0\)

 -~N(\(μ_1-μ_2,{σ_1^2\over n_1}+{σ_2^2\over n_2}\))

两个样本均值差的标准差为

\(se=\sqrt{{σ_1^2\over n_1}+{σ_2^2\over n_2}}=\sqrt{{8^2\over 50}+{8^2\over 50}}=1.6\)

在\(H_0\)为真下，它的正态分布如下

这是一个双尾检验，由上图可知，一旦两个班的均值差落在了-1.96*se=-3.14的左边或者1.96*se=3.14的右边(-1.96和1.96是标准正态分布中的Z值，乘以了se就可以放大到非标准正态分布中来；如果是标准正态分布下，se=1)，就会拒绝\(H_0\)。

如果两班成绩均值差的真值为4分，基于上，问能够检测到4分差异的统计功效是多少？

这里就相当于把\(H_A:μ_1-μ_2≠0\)改成了\(H_A:μ_1-μ_2=4\)

上图中黑线的正态分布就是\(μ_1-μ_2=0\)时的正态分布，绿线的正态分布就是\(μ_1-μ_2=4\)时的正态分布。统计功效就是在3.14右边，绿线正态分布的面积。

Power=P(拒绝\(H_0\) ｜ \(H_A\)为真)，它代表着当两个班的均值差落在了3.14以后的时候，我们会拒绝\(H_0\)，同时要保证\(H_A\)为真时的概率。由于绿线的正态分布与黑线的正态分布的标准差是一样的，都是1.6，则可以计算出绿色区域的面积Power=0.70。

此时我们来看一下α与β之间的关系，α表示在3.14以右，黑线下的面积；而Power就是3.14以右，绿线下的面积。如果我们把3.14这根红色虚线向左移，它意味着增大了α的面积，即增大了犯一类错误的概率，同时也增大了Power的面积，也就是意味着β的减小，也就是犯二类错误的概率在减小。由此我们可以看出如果增大了α，就变相的减小了β；如果我们把3.14这根红色虚线向右移，它意味着减小了α的面积，即减小了犯一类错误的概率，同时也减小了Power的面积，也就是意味着β的增大，也就是犯二类错误的概率在增加。

如果两班成绩均值差的真值为6分，基于上，问能够检测到6分差异的统计功效是多少？

由上图，我们可以很容易求出Power=P(拒绝\(H_0\) ｜ \(H_A\)为真)=0.96。在这里跟上一幅图的区别是绿线的对称轴从4挪到了6，\(H_0\)的参数值和真值之间的距离是会影响Power的，当真值为4的时候，Power=0.70，当真值为6的时候，即拉大了\(H_0\)参数值与真值的距离，Power=0.96。

• 样本容量与统计功效

问样本容量多大可以使得能够检测到4分差异的统计功效是80%？

上图中，红色的虚线依然是0.025黑线下的面积，它距离黑线均值0的距离为1.96*se，它距离绿线均值4的距离为0.84*se，则有

1.96*se+0.84*se=4

se=\(4\over 2.8\)=1.43

\(se=\sqrt{{σ_1^2\over n_1}+{σ_2^2\over n_2}}=\sqrt{{8^2\over n}+{8^2\over n}}=1.43\)

n=63

在计算样本容量时，我们也固定了一些参数，α=0.05，Power=0.80，均值差的真值4，这样我们就可以求出样本容量n了。

1. 当α和样本容量固定时，真值与零假设参数值的间距越大，统计功效越大；
2. 当α和统计功效固定时，真值与零假设参数值的间距越大，所需样本容量越小。

在我们收集数据之前，需要先计算一下样本容量的大小，为了计算这个样本容量，我们通常要设定α，一般设定成0.05，然后需要设定统计功效，一般设定的下限为80%，然后查阅资料来获取真值与零假设之间的间距，然后就可以计算样本容量了。

统计显著性与实际显著性

一班数学成绩的平均分为90，二班数学成绩的平均分为89.5，两班数学成绩分布(总体)的标准差为8，单尾检验

我们假设两个班的样本容量都为n=20，则Z分数为

\(Z={90-89.5\over \sqrt{{8^2\over 20}+{8^2\over 20}}}=0.198\)

p-value=0.42>0.05接受\(H_0\)，我们认为这两个班成绩的均值无差异

我们假设两个班的样本容量都为n=2000，则Z分数为

\(Z={90-89.5\over \sqrt{{8^2\over 2000}+{8^2\over 2000}}}=1.976\)

p-value=0.024<0.05拒绝\(H_0\)，我们认为这两个班成绩的均值有差异

根据以上，我们可以看出样本容量不同的话，我们计算出来的Z值、p值都不同，我们做出来的决策也不同。点估计与零值的真实差异在大样本时更容易检测到(统计显著)。

• 效应量

1. 效应量是量化效应(现象)强度/大小的指标。
2. 效应量有多种测量方式，效应量的绝对值越大表示效应越强，现象越明显。
3. Cohen's D：一种常用的效应量测量，用于量化两个均值差的效应的大小。

\(d={mean(X_1)-mean(X_2)\over S}\)

这里的S为联合标准差

\(S=\sqrt{(n_1-1)S_1^2+(n_2-1)S_2^2\over n_1+n_2-2}\)

回到刚才的数学成绩的示例中，计算d值为

\(d={mean(X_1)-mean(X_2)\over S}={90-89.5\over 8}=0.063<0.2\)

从效应量的角度来看，这是一个效应很小的量，这说明统计显著(statistical significance)≠实际显著(practical significance)，我们在分析问题的时候不仅要关注统计显著，也要关注效应量，两者放到一起考虑才能更好的了解到统计显著性结果有没有实际意义。

Z检验，t检验，卡方检验，F检验的前提条件

• 总体均值

正态总体

1. 组内独立性：样本是通过随机抽样采集的，是随机分配到两个组的。
2. 组间独立性：两组互相独立(非配对)。
3. 如果是无放回抽样，则样本容量不能超过总体容量的10%。

非正态总体

1. 组内和组间独立性。
2. 如果是无放回抽样，则样本容量不能超过总体容量的10%。
3. 样本容量≥30；总体的分布越偏，需要的样本容量越大。

无论是正态总体还是非正态总体，对于总体均值都可以使用Z检验和t检验来进行假设检验，因为中心极限定理告诉我们无论总体分布是正态的还是非正态的，它的样本均值的分布会接近于正态分布。

• 总体方差

正态总体

1. 同总体均值检验的前提条件。可以使用卡方/F检验。

非正态总体

1. 即使满足所有总体均值检验的前提条件，也不能使用卡方/F检验；需要使用非参数方法。

在大部分情况下，我们都关心的是总体均值的情况，对总体方差的情况关注的不是特别大。