---
title: "参数估计-区间估计"
output: 
  html_document:
    number_sections: TRUE
    toc: TRUE
    toc_float: TRUE
    code_folding: show
    fig_width: 7
    fig_height: 5
    theme: united
    highlight: tango
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r,echo=TRUE,message=FALSE,warning =FALSE}

```
<font size=4 face="微软雅黑">

#区间估计

***
&emsp;&emsp;我们之前讲了点估计方法，但是由于点估计只是用一个点去估计未知参数，故在多数情况下会有误差。于是，人们在实际生活中多用区间估计。区间估计是介于估计与检验之间的内容，且区间估计与检验紧密相连。

&emsp;&emsp;先给出区间估计的定义：设$X_1,X_2,...,X_n$为来自参数分布族$\mathcal{F}=\left\{f(x,\theta):\theta\in\Theta\right\}$的样本，$\theta$为一维未知参数。如果$\hat\theta_L(X),\hat\theta_U(X)$为两个统计量,且$\hat\theta_L(X)\leq\hat\theta_U(X)$,则称随机区间$[\hat\theta_L(X),\hat\theta_U(X)]$为$\theta$的一个区间估计。

&emsp;&emsp;从定义中可以看出，给出一个未知参数的区间估计并不难，由于只要给出两个有大小关系的统计量即可。但是，既然是估计，就应有一个好坏的衡量标准，当参数真值为$\theta$时，我们自然希望随机区间$[\hat\theta_L(X),\hat\theta_U(X)]$包含$\theta$的概率$P_{\theta}\left\{\hat\theta_L(X)\leq\theta\leq\hat\theta_U(X)\right\}$越大越好。这个概率就称为置信水平或置信度。由于这个置信水平依赖于参数真值，故我们自然希望对于参数空间中的每一个$\theta$，其置信水平都很大。

##置信区间
<font size=4 face="微软雅黑">

&emsp;&emsp;我们在概率统计中所说的区间估计中的区间指的是置信区间。置信区间的定义是：设$[\hat\theta_L(X),\hat\theta_U(X)]$是参数$\theta$的一个区间估计，如果对给定的$\alpha\in(0,1)$,有$$P_{\theta}\left\{\hat\theta_L(X)\leq\theta\leq\hat\theta_U(X)\right\}\ge1-\alpha,\forall\theta\in\Theta$$

&emsp;&emsp;则称$[\hat\theta_L(X),\hat\theta_U(X)]$为$\theta$的置信水平为$1-\alpha$的置信区间，$\hat\theta_L(X)$和$\hat\theta_U(X)$分别叫置信水平为$1-\alpha$的置信区间的置信下限和置信上限。

&emsp;&emsp;从置信区间的定义可以看出，它可以保证真值落在此区间的概率至少为$1-\alpha$。

##枢轴量法

<font size=4 face="微软雅黑">
&emsp;&emsp;求取参数的置信区间的方法有很多，但其中最常用的是枢轴量法。枢轴量法的基本思想是：在求取某参数的置信区间时，一般先给一个好的点估计，之后再通过这个点估计求感兴趣参数的置信区间。

&emsp;&emsp;用枢轴量法求置信区间分为四步：

* 找一个与待估参数$g(\theta)$无关的统计量$T$（一般是它的一个很好的点估计）；
* 设法找出关于$T$和$g(\theta)$的某函数$S(T,g(\theta))$，使得$S(T,g(\theta))$的分布$F(x)$与$\theta$无关,$S$就称为是枢轴量；
* 适当地选取两个常数$c,d$,使对给定的$\alpha\in(0,1)$，有$P_{\theta}\left\{c\leq S(T,g(\theta))\leq d\right\}=1-\alpha$；
* 如果能把不等式$c\leq S(T,g(\theta))\leq d$等价地改写成$\hat\theta_L(X)\leq g(\theta)\leq\hat\theta_U(X)$，其中，置信上下限只与$c,d,T$有关，而与$\theta$无关，则
$[\hat\theta_L(X),\hat\theta_U(X)]$就是$g(\theta)$的置信水平为$1-\alpha$的置信区间。

#一个正态总体的情况
<font size=4 face="微软雅黑">

***
&emsp;&emsp;用枢轴量法求一个正态总体各参数的置信区间，假设正态总体$X\sim N(\mu,\sigma^2)$,$X_1,X_2,...,X_n$为来自总体$X$的一个样本，$1-\alpha$为置信信度，$\bar{X}$为样本均值，$S^2$为样本方差。
<font size=4 face="微软雅黑">

##均值$\mu$的区间估计
<font size=4 face="微软雅黑">

###$\sigma^2$已知
<font size=4 face="微软雅黑">

&emsp;&emsp;由于$$\bar{X}\sim N(\mu,\sigma_2/n)$$&emsp;&emsp;所以将$\bar{X}$标准化，$$\frac{\bar{X}-\mu}{\sigma/\sqrt{n}}\sim N(0,1)$$&emsp;&emsp;因此有$$P\left\{|\frac{\bar{X}-\mu}{\sigma/\sqrt{n}}|\leq z_{\alpha/2}\right\}=1-\alpha$$其中$z_{\alpha/2}$是正态分布的上$\alpha/2$分位点。

&emsp;&emsp;关于均值$\mu$的置信度为$1-\alpha$的置信区间为$$[\bar{X}-\frac{\sigma}{\sqrt{n}}z_{\alpha/2},\bar{X}+\frac{\sigma}{\sqrt{n}}z_{\alpha/2}]$$
<font size=4 face="微软雅黑">

&emsp;&emsp;在R中，基本统计和数据分析程序包BSDA中已经提供了函数`z.test()`，它可以对基于正态分布的单样本和双样本进行假设检验和区间估计，使用方法如下
```
z.test(x,y=NULL,alternative="two.sided",mu=0,sigma.x=NULL,sigma.y=NULL,conf.level=0.95)
```
&emsp;&emsp;其中，x和y为数值向量，默认y=null，即进行单样本的假设检验；alternative用于指定所求置信区间的类型，默认为two.sided，表示求双尾的值ixnqujian，若为less则求置信上限，为greater求置信下限;mu表示均值，它仅在假设检验中起作用，默认为0；sigma.x和sigma.y分别指定两个样本总体的标准差；conf.level指定区间估计时的置信水平。

&emsp;&emsp;举例说明，首先从均值为10，标准差为2的总体中抽取20个样本，因此这是一个方差已知的正态分布样本。计算置信水平为95%时x的置信区间，调用函数`z.test`
```{r}
set.seed(111)#设定随机种子为确保每次结果一样
x=rnorm(20,10,2)
#install.packages("BSDA")
library(BSDA)
z.test(x,sigma.x=2)$conf.int
```

&emsp;&emsp;结果显示，该样本的置信区间为[8.42,10.17]

###$\sigma^2$未知
<font size=4 face="微软雅黑">

&emsp;&emsp;由于$\sigma^2$未知，所以用样本方差$S^2$代替总体方差$\sigma^2$则$$T=\frac{\bar{X}-\mu}{S/\sqrt{n}}\sim t(n-1)$$有$$P\left\{|\frac{\bar{X}-\mu}{S/\sqrt{n}}|\leq_{\alpha/2}\right\}=1-\alpha$$其中$t_{\alpha/2}(n-1)$表示自由度为$n-1$的$t$分布上的$\alpha/2$分位点。得出关于均值$\mu$置信度为$1-\alpha$的双侧置信区间为$$[\bar{X}-\frac{S}{\sqrt{n}}t_{\alpha/2}(n-1),\bar{X}+\frac{S}{\sqrt{n}}t_{\alpha/2}(n-1)]$$

&emsp;&emsp;在R中实现方差未知的区间估计非常容易，调用t检验的函数就可以直接求得置信区间。函数调用方法如下

```
t.test(x,y=NULL,alternative=c("two.sided","less","greater"),mu=0,paired=FALSE,var.equal=FALSE,conf.level=0.95,...)
```
&emsp;&emsp;其中，x为数据样本；若x和y同时输入，则做爽样本t检验；alternative用于指定所求置信区间的类型，默认为two.sided，表示求双尾的置信区间，若为less则求置信上限，为greater求置信下限；mu表示均值，仅在假设检验中起作用，默认为0。使用上例中的x，假设总体方差未知。
```{r}
t.test(x)$conf.int
```

&emsp;&emsp;总体方差未知时得到x的置信区间为[8.19,10.47]，比方差已知的区间要宽，这个结果十分合理，因为方差未知时我们可用的信息更少，估计结果相对不那么准确。

<font size=4 face="微软雅黑">

##方差$\sigma^2$的区间估计
<font size=4 face="微软雅黑">

###均值$\mu$已知
<font size=4 face="微软雅黑">

&emsp;&emsp;$\sigma^2$的极大似然估计为$$\hat{\sigma}^2=\frac{1}{n}\sum\limits_{i=1}^n(X_i-\mu)^2$$&emsp;&emsp;根据此式可以导出$\sigma^2$的置信区间。
<font size=4 face="微软雅黑">

&emsp;&emsp;$\chi^2(n)$的定义$$\chi^2(n)=\sum\limits_{1=1}^nX_i^2$$其中，$X_i$服从标准正态分布。根据$\chi^2$分布的定义，可以得到$$\frac{n\hat{\sigma}^2}{\sigma^2}=\sum\limits_{i=1}^n(X_i-\mu)^2/\sigma^2\sim \chi^2(n)$$&emsp;&emsp;所以$$P\left\{\chi_{1-\alpha/2}^2(n)\leq\frac{n\hat{\sigma}^2}{\sigma^2}\leq\chi_{\alpha/2}^2(n)\right\}=1-\alpha$$,其中$\chi_{1-\alpha/2}^2(n)$和$\chi_{\alpha/2}^2(n)$分别表示自由度为$n$的$\chi^2$分布上的$1-\alpha/2$和$\alpha/2$分位点。。得到$\sigma^2$的置信度为$1-\alpha$的双侧置信区间为$$[\frac{n\hat{\sigma}^2}{\chi_{\alpha/2}^2(n)},\frac{n\hat{\sigma}^2}{\chi_{1-\alpha/2}^2(n)}]$$
<font size=4 face="微软雅黑">


###均值$\mu$未知
<font size=4 face="微软雅黑">

&emsp;&emsp;取$\sigma^2$的估计量$$S^2=\frac{1}{n-1}\sum\limits_{i=1}^n(X_i-\bar{X})^2$$&emsp;&emsp;并且有$$\frac{(n-1)S^2}{\sigma^2}\sim \chi^2(n-1)$$&emsp;&emsp;因此有$$P\left\{\chi_{1-\alpha/2}^2(n-1)\leq\frac{(n-1)S^2}{\sigma^2}\leq\chi_{\alpha/2}^2(n-1)\right\}=1-\alpha$$,其中$\chi_{1-\alpha/2}^2(n-1)$和$\chi_{\alpha/2}^2(n-1)$分别表示自由度为$n-1$的$\chi^2$分布上的$1-\alpha/2$和$\alpha/2$分位点。由此得到$\sigma^2$的置信度为$1-\alpha$的双侧置信区间为$$[\frac{(n-1)S^2}{\chi_{\alpha/2}^2(n-1)},\frac{(n-1)S^2}{\chi_{1-\alpha/2}^2(n-1)}]$$
<font size=4 face="微软雅黑">

&emsp;&emsp;在R中没有直接计算方差的置信区间的函数，我们可以把上面两种情况写在一个函数里，通过一个if语句进行判断，只要是方差的区间估计，都调用这个函数即可。在R中写函数时，参数可以事先设定一个初值，例如设$mu=Inf$，代表均值未知的情况，调用函数时如果没有特殊说明mu的值，将按照均值未知的方法计算，如果均值已知，在调用函数时对mu重新赋值即可。
```{r}
var.conf.int=function(x,mu=Inf,alpha){
  n=length(x)
  if(mu<Inf){
    s2=sum((x-mu)^2)/n
    df=n
  }
  else{
    s2=var(x)
    df=n-1
  }
  c(df*s2/qchisq(1-alpha/2,df),df*s2/qchisq(alpha/2,df))#qchisq是卡方分布分位数函数
}
```

&emsp;&emsp;在实际分析中，均值$\mu$基本都是未知的情形，例如计算样本X的方差的置信区间时，若均值未知时，调用函数时就不需要指定mu的值了。
```{r}
var.conf.int(x,alpha=0.05)
```

&emsp;&emsp;计算得到总体方差的置信区间是[1.96,7.21]。

#两个正态总体的情况
<font size=4 face="微软雅黑">

***

&emsp;&emsp;单个总体的区间估计针对总体的均值和方差，当我们有两个总体时，就要对它们的参数进行对比了，这在实际应用中很常见。比如，我们要考察一项新技术对提高产品质量是否有成效，则把新技术实施前后的产品质量指标堪称两个正态总体，这时我们所考察的问题，就转化为检验这两个正态总体的均值比较的问题。

##均值差$\mu_1-\mu_2$的区间情况
<font size=4 face="微软雅黑">
&emsp;&emsp;根据两样本总体方差是否已知、是否相等，这里要分三种情况讨论。

###当两总体的方差$\sigma_1^2$,$\sigma_1^2$已知
<font size=4 face="微软雅黑">

&emsp;&emsp;由正态分布的性质有$$\bar{X}-\bar{Y}\sim N(\mu_1-\mu_2,\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2})$$&emsp;&emsp;与单个总体均值区间估计的推导相类似，得到$\mu_1-\mu_2$的置信为$1-\alpha$的双侧置信区间为$$[\bar{X}-\bar{Y}-z_{\alpha/2}\sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}},\bar{X}-\bar{Y}+z_{\alpha/2}\sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}}]$$
<font size=4 face="微软雅黑">

&emsp;&emsp;前面介绍的Z检验函数可以在两总体方差已知的情况下，计算两总体均值差的置信区间，分别用参数sigma.x和sigma.y来说明已知的标准差数值即可。

###当两总体的方差$\sigma_1^2$,$\sigma_1^2$未知但相等
<font size=4 face="微软雅黑">

&emsp;&emsp;经计算可以得到$$T=\frac{\bar{X}-\bar{Y}-(\mu_1-\mu_2)}{S_w\sqrt{1/n_1+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}}$$&emsp;&emsp;推导得到$\mu_1-\mu_2$的置信度为$1-\alpha$的置信区间为$$[\bar{X}-\bar{Y}-t_{\alpha/2}(n_1+n_2-2)S_w\sqrt{1/n_1+1/n_2},\bar{X}-\bar{Y}+t_{\alpha/2}(n_1+n_2-2)S_w\sqrt{1/n_1+1/n_2}]$$

&emsp;&emsp;正如计算单正态总体均值的置信区间，R中的函数$t.test()$还可以用来求两总体均值差的置信区间，由于总体方差相等，需要将其中的参数var.equal设为TRUE。
<font size=4 face="微软雅黑">


###当两总体的方差$\sigma_1^2$,$\sigma_1^2$未知，且不相同
<font size=4 face="微软雅黑">

&emsp;&emsp;可以得到$$T=\frac{\bar{X}-\bar{Y}-(\mu_1-\mu_2)}{\sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}}}\sim t(v)$$其中$$v=(\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2})^2/(\frac{(\sigma_1^2)^2}{n_1^2(n_1-1)}+\frac{(\sigma_2^2)^2}{n_2^2(n_2-1)})$$&emsp;&emsp;但是由于$\sigma_1^2$,$\sigma_2^2$未知,所以可用样本方差$S_1^2$,$S_2^2$来近似，因此$$\hat{v}=(\frac{S_1^2}{n_1}+\frac{S_2^2}{n_2})^2/(\frac{(S_1^2)^2}{n_1^2(n_1-1)}+\frac{(S_2^2)^2}{n_2^2(n_2-1)})$$&emsp;&emsp;所以可以近似的认为$$T\sim t(\hat{v})$$&emsp;&emsp;由此可以得到$\mu_11-\mu_2$的置信度为$1-\alpha$的置信区间为$$[\bar{X}-\bar{Y}-t_{\alpha/2}(\hat{v})\sqrt{\frac{S_1^2}{n_1}+\frac{S_2^2}{n_2}},\bar{X}-\bar{Y}+t_{\alpha/2}(\hat{v})\sqrt{\frac{S_1^2}{n_1}+\frac{S_2^2}{n_2}}]$$
<font size=4 face="微软雅黑">

&emsp;&emsp;这种情况比较复杂，R中没有直接的函数可用，需要手动写出一个函数
```{r}
twosample.ci2=function(x,y,alpha){
  n1=length(x)#计算x长度
  n2=length(y)
  xbar=mean(x)-mean(y)
  s1=var(x)
  s2=var(y)
  nu=(s1/n1+s2/n2)^2/(s1^2/n1^2/(n1-1)+s2^2/n2^2/(n2-1))
  z=qt(1-alpha/2,nu)*sqrt(s1/n1+s2/n2)#t分布分位数函数
  c(xbar-z,xbar+z)
}
```


##方差比$\sigma_1^2/\sigma_2^2$的区间估计
<font size=4 face="微软雅黑">

&emsp;&emsp;在实际分析中，比较两个总体的方差，采用的是“比”的形式，而不是差值，这与统计量的构造有关。

###$\mu_1$与$\mu_2$已知
<font size=4 face="微软雅黑">

&emsp;&emsp;由于$$\hat{\sigma}_1^2=\frac{1}{n_1}\sum\limits_{i=1}^{n_1}(X_i-\mu_1)^2,\hat{\sigma}_2^2=\frac{1}{n_2}\sum\limits_{i=1}^{n_2}(Y_i-\mu_2)^2$$分别为$\sigma_1^2$和$\sigma_2^2$d的无偏估计，由于$$F=\frac{\hat{\sigma_1}^2/\sigma_1^2}{\hat{\sigma_2}^2/\sigma_2^2}\sim F(n_1,n_2)$$&emsp;&emsp;所以$$P\left\{F_{1-\alpha/2}(n_1,n_2)\leq\frac{\hat{\sigma_1}^2/\sigma_1^2}{\hat{\sigma_2}^2/\sigma_2^2}\leq F_{\alpha/2}(n_1,n_2)\right\}=1-\alpha$$，其中$F_{1-\alpha/2}(n_1,n_2)$$ F_{\alpha/2}(n_1,n_2)$/分别为$F(n_1,n_2)$分布上的$1-\alpha/2$和$\alpha/2$分位点。得到$\sigma_1^2/\sigma_2^2$的置信度为$1-\alpha$的置信区间为$$[\frac{\hat{\sigma}_1^2/\hat{\sigma}_2^2}{F_{\alpha/2}(n_1,n_2)},\frac{\hat{\sigma}_1^2/\hat{\sigma}_2^2}{F_{1-\alpha/2}(n_1,n_2)}]$$
<font size=4 face="微软雅黑">

###$\mu_1$与$\mu_2$未知
<font size=4 face="微软雅黑">

&emsp;&emsp;由于$S_1^2$和$S_2^2$分别为$\sigma_1^2$和$\sigma_2^2$的无偏估计，此时有$$F=\frac{S_1^2/\sigma_1^2}{S_2^2/\sigma_2^2}\sim F(n_1-1,n_2-2)$$&emsp;&emsp;因此$$P\left\{F_{1-\alpha/2}(n_1-1,n_2-1)\leq\frac{S_1^2/\sigma_1^2}{S_2^2/\sigma_2^2}\leq F_{\alpha/2}(n_1-1,n_2-1)\right\}=1-\alpha$$则$\mu_1$与$\mu_2$的置信度为$1-\alpha$的置信区间为$$[\frac{S_1^2/S_2^2}{F_{\alpha/2}(n_1-1,n_2-1)},\frac{S_1^2/S_2^2}{F_{1-\alpha/2}(n_1-1,n_2-1)}]$$

&emsp;&emsp;方差比的区间估计与方差的假设检验密不可分，所以R中的函数`var.test()`可以用来直接计算两正态总体方差比的置信区间，调用格式为

```
var.test(x, y, ratio = 1,
         alternative = c("two.sided", "less", "greater"),
         conf.level = 0.95, ...)
```

&emsp;&emsp;其中的参数设置与上文的t检验函数一致，ratio指定原假设中方差比的值，通过alternative设置单侧区间或双侧。

#参考文献
<font size=4 face="微软雅黑">

***
[1]薛毅. 统计建模与R软件[M]. 清华大学出版社, 2007.

[2]王兆军. 数理统计教程[M]. 高等教育出版社,2014.

<font size=4 face="微软雅黑">