统计学part6---参数估计

195 阅读1分钟

image.png

我们开启统计推断的篇章学习啦!我们先从第一个板块参数估计入手,参数估计又分为点估计和区间估计

点估计

估计值、估计量

设总体X的分布函数的形式已知,但它的一一个或多个参数未知,借助于总体X的一个样本来估计总体未知参数的值的问题称为参数的点估计问题

image.png

总体的形式已知(也就是知道总体符合什么分布形式),但是有一个或多个参数未知,那么我们可以在总体中抽取一定的样本,对样本构造统计量,基于统计量获取观察值就认为该观察值近似于未知参数

θ\theta 是待估参数,X1,X2,...,XnX_1,X_2,...,X_n 是总体X的一个样本,x1,x2,...,xnx_1,x_2,...,x_n 是一个样本值。点估计的问题就是要构造一个 适当的统计量 (估计量) ,用它的观察值作为未知参数的近似值 (估计值)

估计量的评选标准

对于同一参数,使用不同的估计方法求出的估计量可能不相同。采用哪一个估计量好呢?

好的估计量需要满足以下几点:

无偏性

若估计量的数学期望存在,并且该期望等于总体参数,则称为无偏估计

均值vs期望:均值是一个统计量(基于样本构造的函数);期望完全由随机变量的概率分布所确定(“上帝视角");两者常混用

对于某些样本值,由这一估计量得到的估计值相对于真值来说偏大或偏小,但是反复将这一估计量使用多次,就“平均”来说其偏差为零。E(估计值)真值E(估计值)-真值称为系统误差;无偏估计的实际意义就是无系统误差.

上面的内容是我们学习样本和抽样分布时所提及的,说明不论总体服从什么分布,样本均值是总体均值的无偏估计;样本方差是总体方差的无偏估计,公式Var(X) = σ2nVar(\overline X) = \frac{\sigma^2}{n} 刻画的是样本均值抽样分布的离散程度,标准误(standard error of mean; SE)

标准误可以理解为是样本均值的标准差(不是样本的标准差!!!)

样本方差:除以n vs除以(n-1)

在描述统计和样本的学习过程中,我们发现方差有时除以n有时除以n-1,下面我们来解释这个问题

统计学part2---描述统计 - 掘金

统计学part5---样本和抽样分布 - 掘金

我们当时介绍了σ 2=i=1n(XiX)2n\sigma ^2=\frac{\displaystyle\sum_{i=1}^n(X_i-\overline{X})^2}{n} S2=i=1n(XiX)2n1S^2=\frac{\displaystyle\sum_{i=1}^n(X_i-\overline{X})^2}{n-1} ,第一个式子中,如果是基于样本计算的,则与总体方差,有系统偏差;第二个式子里面是总体方差的无偏估计

编程实现无偏性

无偏方差和样本方差的差别终于被除数的不同,所以我们可以复用我们之前实现方差的程序并进行修改得到无偏方差的函数

from pq_stats.descriptive_stats import mean, variance
"""无偏性实现"""
def variance_bias(data):
    """有偏方差"""
    n = len(data)
    if n <= 1:
        return None
    mean_value = mean(data)   # 方差的计算需要使用均值
    return sum((e - mean_value) ** 2 for e in data) / n

我们需要抽取很多样本,样本容量还是变化的,对每个样本实现有偏或者无偏的方差的计算,我们通过编写函数实现

类似的函数在中心极限定理时已经实现过了

from pq_stats.descriptive_stats import mean, variance
import random
import matplotlib.pyplot as plt
"""无偏性实现"""
def variance_bias(data):
    """方差"""
    n = len(data)
    if n <= 1:
        return None
    mean_value = mean(data)   # 方差的计算需要使用均值
    return sum((e - mean_value) ** 2 for e in data) / n

def sample(num_of_samples, sample_size, var):  # var是方差函数
    data = []
    for _ in range(num_of_samples):
        # 我们需要对抽取的样本进行操作并将操作的结果放入data中
        data.append(var([random.uniform(0.0, 1.0) for _ in range(sample_size)]))   # 从0~1均匀分布中抽取sample_size个个体组成一个样本,求样本均值
    return data

if __name__ == '__main__':
    # 计算样本的有偏方差
    data1 = sample(1000,40,variance_bias)
    plt.hist(data1, bins='auto', rwidth=0.8)
    plt.axvline(x = mean(data1),c='black')
    plt.axvline(x = 1/12, c='red')
    plt.show()

为什么是112\frac1{12} ?因为# 对于均匀分布来说,总体方差的计算公式为(ba)212\frac{(b-a)^2}{12} ,ab就是均匀分布的两个参数

图中红色的线为理论值,黑色的线为实验值,可以看出有偏的情况下,两者差距还是比较大的,下面我们研究一下无偏的情况:

if __name__ == '__main__':

    # 计算样本的无偏方差
    data2 = sample(1000, 40, variance)
    plt.hist(data2, bins='auto', rwidth=0.8)
    plt.axvline(x=mean(data2), c='black')
    plt.axvline(x=1 / 12, c='red')
    plt.show()

我们看见红线和黑线之间的差距明显比有偏的情况下小了很多,这也证明了在计算方差的时候除以n1n-1 得到的样本方差才是更贴近总体方差的

有效性

有两个无偏估计θ1\theta_1 θ2\theta_2 ,如果在样本容量n相同的情况下,θ1\theta_1 θ2\theta_2 更密集在真值附近,就认为θ1\theta_1 θ2\theta_2 更理想,由于方差是随机变量取值与其数学期望的偏离程度的测量,所以无偏估计以方差最小者为好

image.png

观察上面的图我们可以看出,无论是在左边的图还是右边的图中,绿色的数据点的均值都在靶心处,两者都是无偏估计,但是显然右边的图的数据点偏差比左边的图的数据点大,所以左边的图效果会更好

无偏性和有效性都是在样本容量n固定的前提下提出的

相合性

我们希望随着样本容量的增大,一个估计量的值稳定于待估参数的真值。满足此条件的估计量为相合估计量。

我们假设总体服从标准正态分布N(0,1);样本容量= 20, 70, 120, ... 10000

"""相合性实现"""
import random
from pq_stats.descriptive_stats import mean, variance
import matplotlib.pyplot as plt

if __name__ == '__main__':
    sample_means = []  # 存储样本均值
    sample_vars = []   # 存储样本方差
    indices = []  # 存储每一 个样本均值和样本方差计算它们时候所对应的样本的样本容量

    for sz in range(20,10001,50):
        indices.append(sz)
        sample = [random.gauss(0.0, 1.0) for _ in range(sz)]# 高斯分布抽取sz个样本
        sample_means.append(mean(sample))
        sample_vars.append(variance(sample))

    plt.plot(indices, sample_means)
    plt.plot(indices,sample_vars)
    plt.show()

我们看见代表样本均值(蓝色)的线围绕着0在浮动,当样本容量较小时浮动幅度较大,而随着样本容量的增加,浮动的幅度不断减小;代表样本方差(橘黄色)的线围绕着1在浮动,同样的当样本容量较小时浮动幅度较大,而随着样本容量的增加,浮动的幅度不断减小。通过图我们可以看出我们的样本均值和样本方差逐渐稳定于总体均值和总体方差

区间估计

对于一个未知量,我们在测量或计算时,不仅要得到近似值还要估计误差,即近似值的精确程度/所求真值所在范围,对于未知参数,我们不仅要得到近似值(点估计),还希望估计出一个范围(区间),并希望知道这个范围包含参数真值的可信程度。这种形式的估计称为区间估计,这样的区间称为置信区间

置信区间

设总体的分布函数含有一个未知参数θ,θΘ\theta,\theta\in\Theta ( Θ\Theta θ\theta 可能的取值范围);对于给定值α(0<α<1)\alpha(0\lt\alpha\lt 1) ,若由来自XX的样本X1,X2,...,XnX_1,X_2,...,X_n 确定的两个统计量θ=θ(X1,X2,...,Xn)\underline{\theta}=\underline{\theta}(X_1,X_2,...,X_n) θ=θ(X1,X2,...,Xn)(θ<θ)\overline{\theta}=\overline{\theta}(X_1,X_2,...,X_n)(\underline{\theta}\lt\overline\theta) ,对于任意θΘ\theta\in\Theta 满足P{θ(X1,X2,...,Xn)<θ<θ(X1,X2,...,Xn)}1αP\{\underline{\theta}(X_1,X_2,...,X_n)\lt\theta\lt\overline{\theta}(X_1,X_2,...,X_n)\}\ge1-\alpha ,则称随机区间(θθ)(\underline{\theta},\overline\theta) θθ的置信水平为1α1-\alpha 的置信区间,其中θ\underline{\theta} 称为置信下限,θ\overline\theta 称为置信上限,1α1-\alpha 称为置信水平

固定样本容量n,若反复抽样多次,每个样本值确定一个区间,每个这样的区间要么包含θθ的真值,要么不包含θθ的真值。

image.png

可以看见图中的样本(蓝色的线)与真值要么有交点,要么相离(如圈出来的样本)

按大数定律,在这么多区间中,包含真值的约占100×(1α)%100\times(1-\alpha)\% ,不包含真值的占100×α%100\times\alpha\%

求未知参数θ的置信区间的步骤

  1. 寻求一个样本X1,X2,...,XnX_1,X_2,...,X_n ,和一个统计量W(X1,X2,...,Xn;θ)W(X_1,X_2,...,X_n;\theta) 使统计量W的分布不依赖于θθ和其他未知参数统计量W的构造,通常可以从θθ的点估计着手。
  2. 对于给定的置信水平1α1-\alpha ,定出两个常数a和b,使得P(a<W(X1,X2,...,Xn;θ)< b)=1αP(a\lt W(X_1,X_2,...,X_n;\theta)\lt b)=1-\alpha ,若能从中得到θθ的不等式θ(X1,X2,...,Xn)<θ<θ(X1,X2,...,Xn)\underline{\theta}(X_1,X_2,...,X_n)\lt\theta\lt\overline{\theta}(X_1,X_2,...,X_n) (θθ)(\underline{\theta},\overline\theta) θθ的一个置信水平为1α1-\alpha 的置信区间。

设总体XN(μ,σ2)X\sim N(μ,σ^2) , σ2σ^2 已知,μμ为未知,设X1,X2,...,XnX_1, X_2,...,X_n 是来自X的样本,求μμ的置信水平为1α1-\alpha 的置信区间.

第一步:构建W(X1,X2,...,Xn;θ)W(X_1,X_2,...,X_n;\theta) ,W的分布不依赖于θθ和其他未知参数 Xμσn N(0,1)\frac{\overline {X}-\mu}{\frac{\sigma}{\sqrt{n}}}\sim N(0,1)

第二步:P(a<W(X1,X2,...,Xn;θ)< b)=1αP(a\lt W(X_1,X_2,...,X_n;\theta)\lt b)=1-\alpha

image.png

在标准正态分布中,满足P{X>zα}=α,(0<α<1)P\{X>z_\alpha\}=\alpha,(0<\alpha<1) ,我们需要寻找1α1-\alpha 的位置,那么我们找到点α2\frac\alpha2 的位置,和α2-\frac\alpha2 的位置,两者和标准正态分布曲线组成的面积就是1α1-\alpha

image.png

也就是说a和b的值分别为α2-\frac\alpha2 α2\frac\alpha2

P(a<W(X1,X2,...,Xn;θ)< b)=1α     P{Xμσn<Zα2}=1α     P{XσnZα2<μ<X+σnZα2}=1αP(a\lt W(X_1,X_2,...,X_n;\theta)\lt b)=1-\alpha\implies P\{|\frac{\overline {X}-\mu}{\frac{\sigma}{\sqrt{n}}}|\lt Z_{\frac\alpha2}\}=1-\alpha \implies P\{\overline{X}-\frac{\sigma}{\sqrt{n}}Z_{\frac\alpha2}\lt\mu\lt\overline{X}+\frac{\sigma}{\sqrt{n}}Z_{\frac\alpha2}\}=1-\alpha

其实我们已经得到了μ\mu 的置信水平为1α1-\alpha 的置信区间(XσnZα2,X+σnZα2)(\overline{X}-\frac{\sigma}{\sqrt{n}}Z_{\frac\alpha2},\overline{X}+\frac{\sigma}{\sqrt{n}}Z_{\frac\alpha2}) ,其中σnZα2\frac{\sigma}{\sqrt{n}}Z_{\frac\alpha2} 称作误差范围。

image.png

现在得到的区间属于那些包含μμ的区间的可信程度为100×(1α)%100\times(1-\alpha)\% ,或“该区间包含真值”这一陈述的可信程度为100×(1α)%100\times(1-\alpha)\%

置信水平为1α1-\alpha 的置信区间不是唯一的,选择上面的情况讲解是因为对于正态分布这种单峰且对称的情况,对称的区间长度最短,意味着估计的精度最高

置信区间与样本容量

置信区间(XσnZα2,X+σnZα2)(\overline{X}-\frac{\sigma}{\sqrt{n}}Z_{\frac\alpha2},\overline{X}+\frac{\sigma}{\sqrt{n}}Z_{\frac\alpha2}) ,标准误随着样本容量的增加而减小,误差范围随着样本容量的增加而减小

image.png

一个正态总体的情况

方差已知,求均值的置信区间

已知X N(μ,σ2),X1,X2,...,Xn,1α,X,S2,σ2X\sim N(\mu,\sigma^2),X_1,X_2,...,X_n,1-\alpha,\overline{X},S^2,\sigma^2 ,求均值μ\mu 的置信区间X N(μ,σ2n)\overline{X}\sim N(\mu,\frac{\sigma^2}n) ,Xμσn N(0,1)\frac{\overline {X}-\mu}{\frac{\sigma}{\sqrt{n}}}\sim N(0,1) ,P{Xμσn<Zα2}=1αP\{|\frac{\overline {X}-\mu}{\frac{\sigma}{\sqrt{n}}}|\lt Z_{\frac\alpha2}\}=1-\alpha ,(XσnZα2,X+σnZα2)(\overline{X}-\frac{\sigma}{\sqrt{n}}Z_{\frac\alpha2},\overline{X}+\frac{\sigma}{\sqrt{n}}Z_{\frac\alpha2})

例题

歌曲的时长服从正态分布X N(μ,σ2),σ2=1X\sim N(\mu,\sigma^2),\sigma^2=1 ,手机里有100首歌曲,这100首歌曲的平均时长是4分钟,求μμ的置信水平为95%的置信区间

根据题目可知样本均值=4,样本容量= 100,总体标准差=1,α=0.05\alpha= 0.05 ,(XσnZα2,X+σnZα2)(\overline{X}-\frac{\sigma}{\sqrt{n}}Z_{\frac\alpha2},\overline{X}+\frac{\sigma}{\sqrt{n}}Z_{\frac\alpha2}) ,而95%置信水平对应的Z值为1.96 ,所以(X±σnZα2)= (4±1100×1.96)=(3.804,4.196)(\overline{X}\pm\frac{\sigma}{\sqrt{n}}Z_{\frac\alpha2})= (4\pm\frac1{100}\times1.96)=(3.804,4.196)

Z值通过查表可得

(3.804, 4.196)属于那些包含真值的区间的可信程度为95%; (3.804, 4.196)包含真值这一陈述的可信程度为95%

方差未知,求均值的置信区间

已知X N(μ,σ2),X1,X2,...,Xn,1α,X,S2,X\sim N(\mu,\sigma^2),X_1,X_2,...,X_n,1-\alpha,\overline{X},S^2, 未知σ2\sigma^2 ,求均值μ\mu 的置信区间。

由于σ\sigma 未知,因此Xμσn N(0,1)\frac{\overline {X}-\mu}{\frac{\sigma}{\sqrt{n}}}\sim N(0,1) 不能再适用,我们可以通过用样本方差代替总体方差的方法求解,即XμSn t(n1)\frac{\overline {X}-\mu}{\frac{S}{\sqrt{n}}}\sim t(n-1)

注意:样本容量小时,t分布比正态分布略宽,置信区间略大;这反映了基于小样本得到的结论具有更大的不确定性,随着样本容量增大,t分布近似于标准正态分布;大样本也可以使用标准正态分布来计算置信区间

P{XμSn<tα2(n1)}=1αP\{|\frac{\overline {X}-\mu}{\frac{S}{\sqrt{n}}}|\lt t_{\frac\alpha2}(n-1)\}=1-\alpha (X±Sntα2(n1))(\overline{X}\pm\frac{S}{\sqrt{n}}t_{\frac\alpha2}(n-1))

例题

歌曲的时长服从正态分布X N(μ,σ2),σ2X\sim N(\mu,\sigma^2),\sigma^2 未知,手机里有100首歌曲,这100首歌曲的平均时长是4分钟,方差为1.44,求μμ的置信水平为95%的置信区间

根据题目可知样本均值=4,样本容量= 100,样本标准差=1.2,α=0.05\alpha= 0.05 ,95%置信水平,自由度=n-1= 99,对应的t值为1.98,(X±Sntα2(n1))=4± 1.2100×1.98=(3.762,4.238)(\overline{X}\pm\frac{S}{\sqrt{n}}t_{\frac\alpha2}(n-1))=4\pm \frac{1.2}{\sqrt{100}}\times1.98=(3.762,4.238)

(3.762, 4.238)属于那些包含真值的区间的可信程度为95%; (3.762, 4.238)包含真值这一-陈述的可信程度为95%

均值未知,求方差的置信区间

已知X N(μ,σ2),X1,X2,...,Xn,1α,X,S2,X\sim N(\mu,\sigma^2),X_1,X_2,...,X_n,1-\alpha,\overline{X},S^2, 未知μ\mu ,求求方差 σ2\sigma^2 的置信区间。

这里我们使用(n1)S2σ2 χ2(n1)\frac{(n-1)S^2}{\sigma^2}\sim \chi^2(n-1) ,P{χ1α22(n1)<(n1)S2σ2<χα22(n1)}=1αP\{\chi_{1-\frac\alpha2}^2(n-1)\lt|\frac{(n-1)S^2}{\sigma^2}|\lt \chi_{\frac\alpha2}^2(n-1)\}=1-\alpha ((n1)S2χα22(n1),(n1)S2χ1α22(n1))(\frac{(n-1)S^2}{\chi_{\frac\alpha2}^2(n-1)},\frac{(n-1)S^2}{\chi_{1-\frac\alpha2}^2(n-1)})

例题

歌曲的时长服从正态分布X N(μ,σ2),μX\sim N(\mu,\sigma^2),\mu 未知,手机里有100首歌曲,这100首歌曲的方差为1.44,求σ2\sigma^2 的置信水平为95%的置信区间

根据题目可知样本容量= 100,样本方差=1.2,α=0.05\alpha= 0.05 ((n1)S2χα22(n1),(n1)S2χ1α22(n1))=(99×1.44128.42,99×1.4473.36)=(1.11,1.94)(\frac{(n-1)S^2}{\chi_{\frac\alpha2}^2(n-1)},\frac{(n-1)S^2}{\chi_{1-\frac\alpha2}^2(n-1)})=(\frac{99\times1.44}{128.42},\frac{99\times1.44}{73.36})=(1.11,1.94 )

编程求解

根据上面3种情况的介绍,我们可以发现我们的置信区间要不通过均值去求得,要不通过方差求得,那么我们可以定义两个函数对两者进行求解

from descriptive_stats import mean, std, variance
from math import sqrt
from scipy.stats import norm,t,chi2

def mean_ci_est(data, alpha, sigma=None):
    """估计均值的置信区间"""
    n = len(data)  # 样本容量
    sampl_mean = mean(data)

    if sigma is None:
        # 方差未知
        s = std(data)  # 样本方差
        se = s/sqrt(n)
        t_value = abs(t.ppf(alpha/2, n-1))
        return sampl_mean - se * t_value, sampl_mean + se * t_value

    else:
        # 方差已知
        se = sigma/sqrt(n)  # 标准误
        z_value = abs(norm.ppf(alpha/2))   # 传入面积返回左侧对应的值
        return sampl_mean - se * z_value, sampl_mean + se * z_value


def var_ci_est(data, alpha):
    """求方差的置信区间"""
    n=len(data)
    s2 = variance(data)
    chi2_lower_value = chi2.ppf(alpha/2, n-1)
    chi2_upper_value = chi2.ppf(1-alpha / 2, n - 1)
    
    return (n-1)*s2/chi2_upper_value, (n-1)*s2/chi2_lower_value

两个正态总体的情况

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

已知N(μ1,σ12),X1,X2,...,Xn1,X,S12,N(μ2,σ22),Y1,Y2,...,Yn2,Y,S22,1α,σ12,σ22N(\mu_1,\sigma_1^2),X_1,X_2,...,X_{n_1},\overline{X},S_1^2,N(\mu_2,\sigma_2^2),Y_1,Y_2,...,Y_{n_2},\overline{Y},S_2^2,1-\alpha,\sigma_1^2,\sigma_2^2 ,求均值μ1μ2\mu_1-\mu_2 的置信区间,XY N(μ1μ2,σ12n1+σ22n2)\overline{X}-\overline{Y}\sim N(\mu_1-\mu_2,\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}) ,(XY)(μ1μ2)σ12n1+σ22n2 N(0,1)\frac{(\overline {X}-\overline{Y})-(\mu_1-\mu_2)}{\sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}}}\sim N(0,1) ,(XY±σ12n1+σ22n2Zα2)(\overline{X}-\overline{Y}\pm \sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}}Z_{\frac\alpha2})

例题

25左右人群的月收入服从正态分布N(μ1,σ12)N(μ_1,σ_1^2) , 35左右人群的月收入服从正态分布N(μ2,σ22)N(μ_2, σ_2^2), σ1,σ2\sigma_1,\sigma_2 分别为2000和8000;我们记录了30名25岁和40名35岁个体的月收入。这30名25岁个体平均收入为16,000元,这40名35岁个体平均收入为25,000元。求μ1μ2μ_1-μ_2置信水平为95%的置信区间

(XY±σ12n1+σ22n2Zα2)=(1600025000± 2000230+8000240×1.96)=(10316.56,7683.44)(\overline{X}-\overline{Y}\pm \sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}}Z_{\frac\alpha2})=(16000-25000\pm \sqrt{\frac{2000^2}{30}+\frac{8000^2}{40}}\times1.96)=(-10316.56,-7683.44)

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

已知N(μ1,σ12),X1,X2,...,Xn1,X,S12,N(μ2,σ22),Y1,Y2,...,Yn2,Y,S22,1α,N(\mu_1,\sigma_1^2),X_1,X_2,...,X_{n_1},\overline{X},S_1^2,N(\mu_2,\sigma_2^2),Y_1,Y_2,...,Y_{n_2},\overline{Y},S_2^2,1-\alpha, ,未知σ12=σ22=σ2\sigma_1^2=\sigma_2^2=\sigma^2 求均值μ1μ2\mu_1-\mu_2 的置信区间,(XY)(μ1μ2)Sw1n1+1n2 t(n1+n22),Sw=(n11)S12+(n21)S22n1+n22,(XY±Sw1n1+1n2tα2(n1+n22))\frac{(\overline {X}-\overline{Y})-(\mu_1-\mu_2)}{S_w\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}}\sim t(n_1+n_2-2),S_w=\sqrt{\frac{(n_1-1)S_1^2+(n_2-1)S_2^2}{n_1+n_2-2}},(\overline{X}-\overline{Y}\pm S_w\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}t_{\frac\alpha2}(n_1+n_2-2))

例题

25左右人群的月收入服从正态分布N(μ1,σ12)N(μ_1,σ_1^2) , 35左右人群的月收入服从正态分布N(μ2,σ22)N(μ_2, σ_2^2), σ1,σ2\sigma_1,\sigma_2 相等但未知;我们记录了30名25岁和40名35岁个体的月收入。这30名25岁个体平均收入为16,000元,标准差为2500元,这40名35岁个体平均收入为25,000元,标准差为7000元。求μ1μ2μ_1-μ_2置信水平为95%的置信区间

Sw=(n11)S12+(n21)S22n1+n22=(301)×25002+(401)×7000230+402=5546.925S_w=\sqrt{\frac{(n_1-1)S_1^2+(n_2-1)S_2^2}{n_1+n_2-2}}=\sqrt{\frac{(30-1)\times 2500^2+(40-1)\times7000^2}{30+40-2}}=5546.925

(XY±Sw1n1+1n2tα2(n1+n22))=(9000±5546.925×130+140×1.995)=(11672.72,6327.28)(\overline{X}-\overline{Y}\pm S_w\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}t_{\frac\alpha2}(n_1+n_2-2))=(-9000\pm5546.925\times\sqrt{\frac1{30}+\frac1{40}}\times1.995)=(-11672.72, -6327.28)

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

已知N(μ1,σ12),X1,X2,...,Xn1,X,S12,N(μ2,σ22),Y1,Y2,...,Yn2,Y,S22,1α,σ12σ22N(\mu_1,\sigma_1^2),X_1,X_2,...,X_{n_1},\overline{X},S_1^2,N(\mu_2,\sigma_2^2),Y_1,Y_2,...,Y_{n_2},\overline{Y},S_2^2,1-\alpha,\sigma_1^2\not=\sigma_2^2 ,未知σ12,σ22\sigma_1^2,\sigma_2^2 求均值μ1μ2\mu_1-\mu_2 的置信区间,我们构建Student's t统计量

(XY)(μ1μ2)S12n1+S22n2\frac{(\overline {X}-\overline{Y})-(\mu_1-\mu_2)}{\sqrt{\frac{S_1^2}{n_1}+\frac{S_2^2}{n_2}}} ,该统计量的自由度df=(S12n1+S22n2)2(S12n1)2n11+(S22n2)2n21df=\frac{(\frac{S_1^2}{n_1}+\frac{S_2^2}{n_2})^2}{\frac{(\frac{S_1^2}{n_1})^2}{n_1-1}+\frac{(\frac{S_2^2}{n_2})^2}{n_2-1}} ,(XY±S12n1+S22n2tα2)(\overline{X}-\overline{Y}\pm \sqrt{\frac{S_1^2}{n_1}+\frac{S_2^2}{n_2}}t_{\frac\alpha2})

例题

25左右人群的月收入服从正态分布N(μ1,σ12)N(μ_1,σ_1^2) , 35左右人群的月收入服从正态分布N(μ2,σ22)N(μ_2, σ_2^2), σ1,σ2\sigma_1,\sigma_2 不相等且未知;我们记录了30名25岁和40名35岁个体的月收入。这30名25岁个体平均收入为16,000元,标准差为2500元,这40名35岁个体平均收入为25,000元,标准差为7000元。求μ1μ2μ_1-μ_2置信水平为95%的置信区间

df=(S12n1+S22n2)2(S12n1)2n11+(S22n2)2n21=(2500230+7000240)2(2500230)2301+(7000240)2401=51.394df=\frac{(\frac{S_1^2}{n_1}+\frac{S_2^2}{n_2})^2}{\frac{(\frac{S_1^2}{n_1})^2}{n_1-1}+\frac{(\frac{S_2^2}{n_2})^2}{n_2-1}}=\frac{(\frac{2500^2}{30}+\frac{7000^2}{40})^2}{\frac{(\frac{2500^2}{30})^2}{30-1}+\frac{(\frac{7000^2}{40})^2}{40-1}}=51.394 (XY±S12n1+S22n2tα2)=(9000±2500230+7000240×2.007)=(11402.82,6597.18)(\overline{X}-\overline{Y}\pm \sqrt{\frac{S_1^2}{n_1}+\frac{S_2^2}{n_2}}t_{\frac\alpha2})=(-9000\pm\sqrt{\frac{2500^2}{30}+\frac{7000^2}{40}}\times 2.007)=(-11402.82,-6597.18)

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

已知N(μ1,σ12),X1,X2,...,Xn1,X,S12,N(μ2,σ22),Y1,Y2,...,Yn2,Y,S22,1α,N(\mu_1,\sigma_1^2),X_1,X_2,...,X_{n_1},\overline{X},S_1^2,N(\mu_2,\sigma_2^2),Y_1,Y_2,...,Y_{n_2},\overline{Y},S_2^2,1-\alpha, ,未知μ1,μ2\mu_1,\mu_2 ,求方差比σ12σ22\frac{\sigma_1^2}{\sigma_2^2} 的置信区间

S12/S22σ12/σ22F(n11,n21)\frac{{S_1^2}/{S_2^2}}{{\sigma_1^2}/{\sigma_2^2}}\sim F(n_1-1,n_2-1) ,(S12S221Fα2(n11,n21),S12S221F1α2(n11,n21))(\frac{S_1^2}{S_2^2}\frac1{F_{\frac\alpha2}(n_1-1,n_2-1)},\frac{S_1^2}{S_2^2}\frac1{F_{1-\frac\alpha2}(n_1-1,n_2-1)})

例题

25左右人群的月收入服从正态分布N(μ1,σ12)N(μ_1,σ_1^2) , 35左右人群的月收入服从正态分布N(μ2,σ22)N(μ_2, σ_2^2), μ1,μ2\mu_1,\mu_2 未知;我们记录了30名25岁和40名35岁个体的月收入。这30名25岁个体平均收入为16,000元,标准差为2500元,这40名35岁个体平均收入为25,000元,标准差为7000元。求σ12σ22\frac{\sigma_1^2}{\sigma_2^2} 置信水平为95%的置信区间

(S12S221Fα2(n11,n21),S12S221F1α2(n11,n21))=(2500270002×11.962,2500270002×10.492)=(0.065,0.259)(\frac{S_1^2}{S_2^2}\frac1{F_{\frac\alpha2}(n_1-1,n_2-1)},\frac{S_1^2}{S_2^2}\frac1{F_{1-\frac\alpha2}(n_1-1,n_2-1)})=(\frac{2500^2}{7000^2}\times\frac1{1.962},\frac{2500^2}{7000^2}\times\frac1{0.492})=(0.065,0.259)

编程实现

根据上面的介绍我们可以知道,前两种情况都是方差未知的情况,后两种情况是方差已知的情况,我们先从方差已知的情况入手

def mean_diff_ci_t_est(data1, data2, alpha, equal = True):  # equal代表方差是否相等
    """两个总体方差未知,求均值差的置信区间"""
    n1 = len(data1)
    n2 = len(data2)
    mean_diff = mean(data1) - mean(data2)  # 均值差

    sample1_var = variance(data1)  # 方差
    sample2_var = variance(data2)

    if equal:
        sw = sqrt(((n1 - 1) * sample1_var + (n2 - 1) * sample2_var) / (n1 + n2 - 2))
        t_value = abs(t.ppf(alpha / 2, n1 + n2 - 2))
        return mean_diff - sw * sqrt(1 / n1 + 1 / n2) * t_value, mean_diff + sw * sqrt(1 / n1 + 1 / n2) * t_value

    else:
        df_numerator = (sample1_var / n1 + sample2_var / n2) ** 2
        df_denominator = (sample1_var / n1) ** 2 / (n1 - 1) + (sample2_var / n2) ** 2 / (n2 - 1)
        df = df_numerator / df_denominator
        t_valus = abs(t.ppf(alpha / 2, df))
        return mean_diff - sqrt(sample1_var / n1 + sample2_var / n2) * t_valus, mean_diff + sqrt(
            sample1_var / n1 + sample2_var / n2) * t_valus

def main_diff_ci_z_est(data1, data2, alpha, sigma1, sigma2):
    """两个总体方差已知,求均值差的置信区间"""

    n1 = len(data1)
    n2 = len(data2)

    mean_diff = mean(data1) - mean(data2)
    z_value = abs(norm.ppf(alpha / 2))

    return mean_diff - sqrt(sigma1 ** 2 / n1 + sigma2 ** 2 / n2) * z_value, mean_diff + sqrt(
        sigma1 ** 2 / n1 + sigma2 ** 2 / n2) * z_value

def var_ratio_ci_est(data1, data2, alpha):
    """两个总体均值未知,求方差比的置信区间"""
    n1 = len(data1)
    n2 = len(data2)
    f_lower_value = f.ppf(alpha / 2, n1 - 1, n2 - 1)
    f_upper_value = f.ppf(1 - alpha / 2, n1 - 1, n2 - 1)
    var_ratio = variance(data1) / variance(data2)
    return var_ratio / f_upper_value, var_ratio / f_lower_value

单侧置信区间

在某些实际问题中,我们只关心“上限"或者“下限”

电池/灯泡的平均寿命;有害物质的平均含量

P{θ(X1,X2,...,Xn)<θ}1α (θ,+)P\{\underline{\theta}(X_1,X_2,...,X_n)\lt\theta\}\ge1-\alpha \to(\underline{\theta},+\infty) 单侧置信下限

P{θ<θ(X1,X2,...,Xn)}1α(,θ)P\{\theta\lt\overline{\theta}(X_1,X_2,...,X_n)\}\ge1-\alpha\to(-\infty,\overline{\theta}) 单侧置信上限

P为置信水平为1α1-\alpha 的单侧置信区间

总结

image.png