python 实现方差检验

发布于:2024-04-20 ⋅ 阅读:(32) ⋅ 点赞:(0)

前面介绍的 z z z检验和 t t t检验主要用于对总体的均值进行假设检验,下面主要介绍对总体方差进行假设检验的方法。

一. 单个总体方差检验

设总体 X ∼ N ( μ , δ 2 ) X\sim N(\mu, \delta^2) XN(μ,δ2) μ , δ 2 \mu,\delta^2 μ,δ2均未知, X 1 , X 2 . . . , X n X_1,X_2...,X_n X1,X2...,Xn是来自 X X X的样本,则
( n − 1 ) S 2 δ 2 ∼ χ 2 ( n − 1 ) \frac{(n- 1)S^2}{\delta^2} \sim \chi^2(n - 1) δ2(n1)S2χ2(n1)

其中S为样本标准差, δ \delta δ为总体标准差,n为样本容量

根据以上的统计量,可以对单个总体的方差 δ \delta δ进行假设检验。

原假设 H 0 H_0 H0 备则假设 H 1 H_1 H1 拒绝域
δ 2 = δ 0 2 \delta^2 = \delta_0^2 δ2=δ02 δ 2 ≠ δ 0 2 \delta^2 \neq \delta_0^2 δ2=δ02 χ 2 > χ α / 2 2 ( n − 1 ) \chi^2 > \chi_{\alpha/2}^2(n - 1) χ2>χα/22(n1) χ 2 < χ 1 − α / 2 2 ( n − 1 ) \chi^2 < \chi_{1- \alpha/2}^2(n - 1) χ2<χ1α/22(n1)
δ 2 ≤ δ 0 2 \delta^2 \leq \delta_0^2 δ2δ02 δ 2 > δ 0 2 \delta^2 > \delta_0^2 δ2>δ02 χ 2 > χ α 2 ( n − 1 ) \chi^2 > \chi_{\alpha}^2(n - 1) χ2>χα2(n1)
δ 2 ≥ δ 0 2 \delta^2 \geq \delta_0^2 δ2δ02 δ 2 < δ 0 2 \delta^2 < \delta_0^2 δ2<δ02 χ 2 < χ 1 − α 2 ( n − 1 ) \chi^2 < \chi_{1 - \alpha}^2(n - 1) χ2<χ1α2(n1)

某厂生产某种型号的电池,其寿命长期以来服从方差 δ 2 = 5000 \delta^2 = 5000 δ2=5000的正太分布,现有一批这种电池,从它的生产情况看,寿命的波动性有所改变。现随机抽26只电池,测出其寿命的样本方差 S 2 = 9200 S^2 = 9200 S2=9200,请问根据这一数据能否推断这批电池寿命的波动性较以前有显著的改变?

1. 双侧检验

α = 0.05 \alpha = 0.05 α=0.05的水平下做出如下假设:
H 0 : δ 2 = 5000    H 1 : δ 2 ≠ 5000 H_0: \delta^2 = 5000 \space\space H_1: \delta^2 \neq 5000 H0:δ2=5000  H1:δ2=5000

从备则假设可知,该假设检验是一个双边检验。

代码实现:

from scipy.stats import chi2

if __name__ == '__main__':

    alpha = 0.05
    n = 26

    # 总体方差
    var = 5000
    # 样本方差
    sample_var = 10000
    # 卡方统计量
    chi2stats = (n - 1) * var / sample_var

    rv = chi2(df=n - 1)
    chi2left = rv.ppf(alpha / 2)
    chi2right = rv.ppf(1 - alpha / 2)

    print("chi2left: {}, chi2right: {}".format(round(chi2left, 3), round(chi2right, 3)))
    print("chi2stats: ", round(chi2stats, 3))

    # 计算p值
    if chi2stats > chi2right:
        pval = 2 * (1 - rv.cdf(chi2stats))
    else:
        pval = 2 * (rv.cdf(chi2stats))

    if chi2stats > chi2right or chi2stats < chi2left:
        print("reject null hypothesis, pval: ", round(pval, 3))
    else:
        print("not reject null hypothesis, pval: ", round(pval, 3))

运行结果:

chi2left: 13.12, chi2right: 40.646
chi2stats:  12.5
reject null hypothesis, pval:  0.036
2. 单侧检验

保持 α = 0.05 \alpha = 0.05 α=0.05不变,修改原假设和备则假设 H 0 : δ 2 ≥ 5000   H 1 : δ 2 < 5000 H_0:\delta^2 \geq 5000 \space H_1: \delta^2 <5000 H0:δ25000 H1:δ2<5000,此时就变成了一个左侧检验。

代码实现:

from scipy.stats import chi2

if __name__ == '__main__':

    alpha = 0.05
    n = 26

    # 总体方差
    var = 5000
    # 样本方差
    sample_var = 10000
    # 卡方统计量
    chi2stats = (n - 1) * var / sample_var

    rv = chi2(df=n - 1)
    chi2left = rv.ppf(alpha)

    print("chi2left: {}".format(round(chi2left, 3)))
    print("chi2stats: ", round(chi2stats, 3))

    # 计算p值
    pval = rv.cdf(chi2stats)

    if chi2stats < chi2left:
        print("reject null hypothesis, pval: ", round(pval, 3))
    else:
        print("not reject null hypothesis, pval: ", round(pval, 3))

运行结果:

chi2left: 14.611
chi2stats:  12.5
reject null hypothesis, pval:  0.018
二. 两个总体方差检验

X 1 , X 2 . . . , X n X_1,X_2...,X_n X1,X2...,Xn是来自总体 N ( μ 1 , δ 1 2 ) N(\mu_1,\delta_1^2) N(μ1,δ12)的样本, Y 1 , Y 2 . . . , Y n Y_1,Y_2...,Y_n Y1,Y2...,Yn是来自总体 N ( μ 2 , δ 2 2 ) N(\mu_2,\delta_2^2) N(μ2,δ22)的样本,且两样本独立,其样本方差分别为 S 1 2 , S 2 2 S_1^2,S_2^2 S12,S22,且 μ 1 , μ 2 , δ 1 2 , δ 2 2 \mu_1,\mu_2,\delta_1^2, \delta_2^2 μ1,μ2,δ12,δ22均未知,则
S 1 2 / S 2 2 δ 1 2 / δ 2 2 ∼ F ( n 1 − 1 , n 2 − 1 ) \frac{S_1^2/S_2^2}{\delta_1^2/\delta_2^2} \sim F(n_1 - 1, n_2 - 1) δ12/δ22S12/S22F(n11,n21)

δ 1 , δ 2 , \delta_1, \delta_2, δ1,δ2,为两个总体的标准差, n 1 , n 2 n_1,n_2 n1,n2为两个样本的容量

根据以上的统计量,可以对两个总体的方差之间的大小关系 δ 1 2 / δ 2 2 \delta_1^2/\delta_2^2 δ12/δ22进行假设检验。

原假设 H 0 H_0 H0 备则假设 H 1 H_1 H1 拒绝域
δ 1 2 = δ 2 2 \delta_1^2 = \delta_2^2 δ12=δ22 δ 1 2 ≠ δ 2 2 \delta_1^2 \neq \delta_2^2 δ12=δ22 F > F α / 2 ( n 1 − 1 , n 2 − 1 ) F > F_{\alpha/2}(n_1 - 1, n_2 - 1) F>Fα/2(n11,n21) F < F 1 − α / 2 ( n 1 − 1 , n 2 − 1 ) F < F_{1 - \alpha/2}(n_1 - 1, n_2 - 1) F<F1α/2(n11,n21)
δ 1 2 ≤ δ 2 2 \delta_1^2 \leq \delta_2^2 δ12δ22 δ 1 2 > δ 2 2 \delta_1^2 > \delta_2^2 δ12>δ22 F > F α ( n 1 − 1 , n 2 − 1 ) F > F_{\alpha}(n_1 - 1, n_2 - 1) F>Fα(n11,n21)
δ 1 2 ≥ δ 2 2 \delta_1^2 \geq \delta_2^2 δ12δ22 δ 1 2 < δ 2 2 \delta_1^2 < \delta_2^2 δ12<δ22 F < F 1 − α ( n 1 − 1 , n 2 − 1 ) F <F_{1 - \alpha}(n_1 - 1, n_2 - 1) F<F1α(n11,n21)

两工厂生产同一份零件且零件尺寸服从正态分布,A生产的16个零件中,均值 μ A \mu_A μA为3.1,样本方差 δ A 2 \delta_A^2 δA2为0.04;B生产的25个零件中,均值 μ B \mu_B μB为3.08,样本方差 δ B 2 \delta_B^2 δB2为0.09。问能否在 α = 0.05 \alpha = 0.05 α=0.05的显著性水平下认为两工厂生产零件的方差相同?

1. 双边检验

α = 0.05 \alpha = 0.05 α=0.05的水平下做出如下假设:
H 0 : δ A 2 = δ B 2    H 1 : δ A 2 ≠ δ B 2 H_0: \delta_A^2 = \delta_B^2 \space\space H_1: \delta_A^2 \neq \delta_B^2 H0:δA2=δB2  H1:δA2=δB2
从备则假设可知,该假设检验是一个双边检验。

代码实现:

from scipy.stats import f

if __name__ == '__main__':

    # 两样本的方差与样本数量
    var1, var2 = 0.04, 0.09
    df1, df2 = 16, 25

    alpha = 0.05
    rv = f(dfn= df1 - 1, dfd= df2 - 1)
    # 计算拒绝域
    fleft = rv.ppf(alpha / 2)
    fright = rv.ppf(1 - alpha / 2)
    # 计算统计量
    fstats = var1 / var2

    print("fleft: {}, fright: {}".format(round(fleft, 3), round(fright, 3)))
    print("fstats: ", round(fstats, 3))
    pval = rv.sf(fstats)

    if fstats > fright:
        pval = 2 * rv.sf(fstats)
    else:
        pval = 2 * rv.cdf(fstats)

    if fstats > fright or fstats < fleft:
        print("reject null hypothesis, pval is ", round(pval, 3))
    else:
        print("not reject null hypothesis, pval is ", round(pval, 3))

运行结果:

fleft: 0.37, fright: 2.437
fstats:  0.444
not reject null hypothesis, pval is  0.107
2. 单边检验

保持 α = 0.05 \alpha = 0.05 α=0.05不变,修改原假设和备则假设 H 0 : δ A 2 ≥ δ B 2   H 1 : δ A 2 < δ B 2 H_0:\delta_A^2 \geq \delta_B^2 \space H_1: \delta_A^2 < \delta_B^2 H0:δA2δB2 H1:δA2<δB2,此时就变成了一个左侧检验。

代码实现:

from scipy.stats import f

if __name__ == '__main__':
    # 两样本的方差与样本数量
    var1, var2 = 0.04, 0.09
    df1, df2 = 16, 25

    alpha = 0.05
    rv = f(dfn=df1 - 1, dfd=df2 - 1)
    # 计算拒绝域
    fleft = rv.ppf(alpha)

    # 计算统计量
    fstats = var1 / var2

    print("fleft: {}".format(round(fleft, 3)))
    print("fstats: ", round(fstats, 3))
    pval = rv.cdf(fstats)


    if fstats < fleft:
        print("reject null hypothesis, pval is ", round(pval, 3))
    else:
        print("not reject null hypothesis, pval is ", round(pval, 3))

运行结果:

fleft: 0.437
fstats:  0.444
not reject null hypothesis, pval is  0.053