R语言学习笔记（四）参数估计

文章目录

写在前面

这次总结一下数理统计中的参数估计，即**点估计(矩估计、极大似然估计)和区间估计(置信区间)**部分的​​R​​语言实现，由于这部分内容没有相应的R语言内置函数，所以需要编程的地方比较多，篇幅也相应地比较长。

• 在计算非线性方程组的根时，采用了自定义函数​​Newtons()​​，运用Newton法进行求根。

# 定义Newton法迭代的函数：计算非线性方程组
Newtons <- function(fun, x, eps=1e-5, it_max=100) {
    index <- 0; k <- 1;
    while (k <= it_max) {
        x1 <- x; obj <- fun(x);
        x <- x - solve(obj$J, obj$f);
        norm <- sqrt((x - x1) %*% (x - x1))
        # 达到精度，跳出循环，index赋值为1表示计算成功
        if (norm < eps) {
            index <- 1; break
        }
        k <- k + 1
    }
    obj <- fun(x);
    list(root=x, it_num=k, index=index, FunVal=obj$f)
}

点估计

极大似然估计

极大似然估计(Maximum Likelihood Estimate, MLE)，最早由统计学家Fisher提出，是一种充分利用总体分布函数信息的估计方式，方法是寻找使似然函数达到最大的参数$\theta$。

• 定义：设总体X的概率密度函数或分布律为$f(x;\,\theta),\,\theta\in\Theta$是未知参数，$X_1,\,X_2,\,\cdots,\,X_n$为来自总体$X$的样本，称
  $L(\theta;\,x)=L(\theta;x_1,\,x_2,\,\cdots,\,x_n)=\prod\limits_{i=1}^nf(x_i;\,\theta)$
  为$\theta$的极大似然函数(likelihood function)。
• 定义：设总体X的概率密度函数或分布律为$f(x;\,\theta),\,\theta\in\Theta$是未知参数，$X_1,\,X_2,\,\cdots,\,X_n$为来自总体$X$的样本，$L(\theta;\,x)$为$\theta$的似然函数, 若$\hat{\theta}=\hat{\theta}(X)=\hat{\theta}(X_1,\,X_2,\,\cdots,\,X_n)$是一个统计量，且满足：
  $L(\hat{\theta}(X);\,X)=\sup\limits_{\theta\in\Theta}L(\theta;\,X)$
  则称$\hat{\theta}$为$\theta$的最大似然估计。

下面介绍几种常见分布的似然函数及其推导。

• 均匀分布
  显然得到$\hat{a}=X_{(1)},\,\hat{b}=X_{(n)}$.
• 指数分布
  服从指数分布的最大似然估计函数为
  $L(\lambda;\,x) =\lambda^n\mathrm{e}^{-\lambda\sum\limits_{i=1}^nx_i}$
  取对数并求导得到
  $\frac{\partial \ln L(\lambda;\,x)}{\partial \lambda} =\left(n\ln\lambda-\lambda\sum\limits_{i=1}^nx_i\right)_{\lambda} =\frac{n}{\lambda}-\sum_{i=1}^n x_i=0$
  即$\lambda=\dfrac{n}{\sum\limits_{i=1}^nx_i}$.
• 正态分布
  正态分布的似然函数为

$L(\mu,\,\sigma^2;\,x)=\prod_{i=1}^nf(x_i;\,\mu,\,\sigma^2)=(2\pi\sigma^2)^{-\frac n2}\exp\left[-\frac1{2\sigma^2}\sum_{i=1}^n(x_i-\mu)^2\right],$

对数似然函数为
$\ln L(\mu,\,\sigma^2;\,x) =-\frac n2\ln(2\pi\sigma^2)-\frac1{2\sigma^2}\sum_{i=1}^n(x_i-\mu)^2,$
令
$\begin{cases} \dfrac{\partial \ln L(\mu,\,\sigma^2;\,x)}{\partial \mu} =\dfrac1{\sigma^2}\sum\limits_{i=1}^n(x_i-\mu)=0, \\ \dfrac{\partial \ln L(\mu,\,\sigma^2;\,x)}{\partial \sigma^2} =-\dfrac{n}{2\sigma^2}+\dfrac1{2\sigma^4}\sum\limits_{i=1}^n(x_i-\mu)^2=0, \end{cases}$
解此似然方程组得到：
$\mu=\dfrac1n\sum\limits_{i=1}^nx_i=\overline{x},\quad \sigma^2=\dfrac1n\sum_{i=1}^n(x_i-\overline{x})^2,$
进一步验证，对于对数似然函数$\ln L(\mu,\,\sigma^2;\,x)$的二阶Hesse矩阵
$\begin{bmatrix} -\dfrac n{\sigma^2} & -\dfrac1{\sigma^4}\sum\limits_{i=1}^n(x_i-\mu)\\ -\dfrac1{\sigma^4}\sum\limits_{i=1}^n(x_i-\mu) & \dfrac n{2\sigma^4}-\dfrac1{\sigma^6}\sum\limits_{i=1}^n(x_i-\mu)^2 \end{bmatrix} = \begin{bmatrix} -\dfrac n{\sigma^2} & 0\\ 0 & -\dfrac{n}{2\sigma^4} \end{bmatrix}$
为负定矩阵，所以$\left(\overline{x},\,\dfrac1n\sum\limits_{i=1}^n(x_i-\overline{x})^2\right)$是$L(\mu,\,\sigma^2;\,x)$的极大值点。故$(\mu,\,\sigma^2)$的最大似然估计为
$\hat{\mu}=\overline{X},\quad\hat{\sigma}^2=\dfrac1n\sum\limits_{i=1}^n(X_i-\overline{X})^2.$

下面分两种情况进行极大似然估计中参数的计算。

可求出解析解

首先采用Newton法实现：

# 定义待求方程
model <- function(e) {
    set.seed(7)
    x <- rnorm(10)
    n <- length(x)
    f <- c(sum(x - e[1]), 
           -n + sum((x - e[1])^2 / e[2]^4))
    J <- matrix(c(-n, 0, -2 * sum(x - e[1]) / e[2]^4,
                -4 * e[2]^(-3) * sum((x - e[1])^2)), 
                nrow=2, byrow=T)
    list(f=f, J=J)
}

# 调用自定义函数`Newtons()`进行求解
Newtons(model, c(0, 1))

## $root
## [1] 0.1039757 1.0962031
## 
## $it_num
## [1] 7
## 
## $index
## [1] 1
## 
## $FunVal
## [1] -3.608225e-16  1.878941e-05

下面介绍一个简单的方法，需要调用​​rootSolve​​​外部包的​​multiroot()​​​函数，求解有$n$个方程、$n$个未知量的非线性方程组。

# 定义待求方程
model <- function(e, x) {
    n <- length(x)
    F1 <- sum(x - e[1])
    F2 <- -n + sum((x - e[1])^2 / e[2]^4)
    c(F1, F2)
}

# 调用函数`multiroot()`进行求解
set.seed(7)
x <- rnorm(10)
# 导入外部包
library(rootSolve)
#  求解
multiroot(f=model, start=c(0, 1), x=x)

## $root
## [1] 0.1039757 1.0962036
## 
## $f.root
## [1] -3.469447e-16  5.412950e-10
## 
## $iter
## [1] 5
## 
## $estim.precis
## [1] 2.706477e-10

不易或无法求出解析解

采用数值解法

以Cauchy分布的最大似然估计为例

• 采用​​uniroot()​​函数

# 参数为1的cauchy分布
set.seed(7)
x <- rcauchy(100,　1)
f <- function(p) sum((x - p) / (1 + (x - p)^2))
out <- uniroot(f, c(0, 5)); out

## $root
## [1] 1.08361
## 
## $f.root
## [1] -0.0001693485
## 
## $iter
## [1] 6
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05

• 采用​​optimize()​​​函数，可以达到与​​uniroot()​​函数一致的结果。

# 生成参数为1的Cauchy分布样本
set.seed(7)
x <- rcauchy(100, 1)
loglike <- function(p) {
    n <- length(x)
    -log(pi) * n - sum(log(1 + (x - p)^2))
}
optimize(loglike, c(0, 5), maximum = T)

## $maximum
## [1] 1.083612
## 
## $objective
## [1] -257.9063

矩估计

使用矩估计进行参数估计的方法称为矩法(method of moments)，由英国统计学家K · Pearson提出，思想是用样本矩去估计总体矩，总体矩与总体的参数有关，从而得到总体参数的估计。

利用矩法估计总体的均值和方差，就等价于用样本的一阶原点矩估计均值，用样本的二阶中心矩估计方差。

下面介绍一些常用分布的矩估计推导。

• 均匀分布
  分为两种情况，第一种只需要求解一阶原点矩，而第二种(一般情况)还需要计算二阶中心矩。

• 情形一（特殊情况）
  $EX=\int_0^\theta x\frac1\theta\mathrm{d}x=\frac\theta2,$
  所以其矩估计为$\hat{\theta}=2\overline{X}=\dfrac2n\sum\limits_{i=1}^nX_i$.
• 情形二（一般情况）

$\begin{aligned} EX&=\int_a^b x\frac1{b-a}\mathrm{d}x=\frac{b+a}2,\\ DX&=\int_a^b x^2\frac1{b-a}\mathrm{d}x-\left(\frac{b+a}2\right)^2=\frac{(b-a)^2}{12}, \end{aligned}$
令
$\begin{cases} \dfrac{b+a}2=\overline{X}\\ \dfrac{(b-a)^2}{12}=\dfrac1n\sum\limits_{i=1}^nX_i^2 \end{cases}$
解得$\hat{a}=\overline{X}-\sqrt{\dfrac3n\sum\limits_{i=1}^nX_i^2},\ \hat{b}=\overline{X}+\sqrt{\dfrac3n\sum\limits_{i=1}^nX_i^2}.$

• 指数分布
  $EX=\int_0^{+\infty}\lambda x\mathrm{e}^{-\lambda x}\mathrm{d}x=\frac1\lambda,$
  因此其矩估计为$\hat{\lambda}=\dfrac{n}{\sum\limits_{i=1}^{n}X_i}$.
• 正态分布
  算总体$X$的一阶、二阶原点矩
  $M_1 =EX=\mu,\quad M_2 =EX^2=\sigma^2+\mu^2$
  以及样本的一阶、二阶原点矩
  $A_1=\overline{X}=\frac1n\sum_{i=1}^nX_i,\quad A_2=\frac1n\sum_{i=1}^nX_i^2.$
  所以得到方程组
  $\begin{cases} \mu=\overline{X}\\ \sigma^2+\mu^2=\dfrac1n\sum\limits_{i=1}^nX_i^2 \end{cases}$
  解上述方程，得均值$\mu$和方差$\sigma^2$的矩估计
  $\begin{aligned}\hat{\mu}&=\overline{X},\\\hat{\sigma}^2&=\dfrac1n\sum\limits_{i=1}^nX_i^2-\overline{X}^2=\dfrac1n\sum\limits_{i=1}^n(X_i-\overline{X})^2=\frac{n-1}nS^2.\end{aligned}$

# 定义待求方程组
moment_fun <- function(p) {
    f <- c(p[1] * p[2] - M1, 
           p[1] * p[2] - p[1] * p[2]^2 - M2)
    J <- matrix(c(p[2], p[1], p[2] - p[2]^2, 
                  p[1] - 2 * p[1] * p[2]), nrow=2, byrow=T)
    list(f=f, J=J)
}


# 主函数
# N=20, p=0.7, 试验次数n=100
# 设置随机数种子，使每次运行得到相同的结果
set.seed(7)
# 生成服从二项分布的随机数作为输入数据
x <- rbinom(100, 20, 0.7);
n <- length(x)
M1 <- mean(x)
M2 <- (n-1) / n * var(x)
# 计算矩估计参数
p <- c(10, 0.5); 
Newtons(moment_fun, p)

## $root
## [1] 20.1441323  0.6875451
## 
## $it_num
## [1] 6
## 
## $index
## [1] 1
## 
## $FunVal
## [1] -1.776357e-15  8.881784e-16

区间估计

这部分的内容比较多，因为涉及到的情况分类多。不过编程不难，直接根据公式与对应的适应情况进行编程即可，主要用到了​​if-else​​条件分支语句。

一个正态总体的置信区间

$\sigma^2$已知时，$\mu$的区间估计

# 编写函数计算置信区间
# sigma默认取值为-1，代表sigma未知的情况
interval_estimate1 <- function(x, sigma=-1, alpha=0.05) { 
   n <- length(x); 
   xb <- mean(x)
   if (sigma >= 0) {
     tmp <- sigma / sqrt(n) * qnorm(1 - alpha / 2); 
     df <- n
     }
   else {
     tmp <- sd(x) / sqrt(n) * qt(1 - alpha / 2, n - 1);
     df <- n - 1   
     }
   list(mean=xb, df=df, a=xb - tmp, b=xb + tmp)
}

# 例题求解
x <- c(14.6, 15.1, 14.9, 14.8, 15.2, 15.1)
interval_estimate1(x, sigma=0.2)

## $mean
## [1] 14.95
## 
## $df
## [1] 6
## 
## $a
## [1] 14.78997
## 
## $b
## [1] 15.11003

t.test(x)

## 
##  One Sample t-test
## 
## data:  x
## t = 162.16, df = 5, p-value = 1.692e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  14.713 15.187
## sample estimates:
## mean of x 
##     14.95

$\sigma^2$未知时，$\mu$的区间估计

interval_estimate1(x)

## $mean
## [1] 14.95
## 
## $df
## [1] 5
## 
## $a
## [1] 14.713
## 
## $b
## [1] 15.187

t.test(x)

## 
##  One Sample t-test
## 
## data:  x
## t = 162.16, df = 5, p-value = 1.692e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  14.713 15.187
## sample estimates:
## mean of x 
##     14.95

方差$\sigma^2$的区间估计

# 编写自定义函数计算置信区间
# 默认mu=Inf，代表mu未知的情况
interval_var1 <- function(x, mu=Inf, alpha=0.05) { 
   n <- length(x) 
   if (mu < Inf) {
     S2 <- sum((x - mu)^2) / n; 
     df <- n   
     }
   else{      
     S2 <- var(x); 
     df <- n-1
     }
   a <- df * S2 / qchisq(1 - alpha / 2, df)
   b <- df * S2 / qchisq(alpha / 2, df)
   list(var=S2, df=df, a=a, b=b)
}

# 例题求解
x <- c(10.1, 10, 9.8, 10.5, 9.7, 10.1, 9.9, 10.2, 10.3, 9.9)

# mu已知
interval_var1(x, mu=10)

## $var
## [1] 0.055
## 
## $df
## [1] 10
## 
## $a
## [1] 0.0268513
## 
## $b
## [1] 0.1693885

# mu未知
interval_var1(x)

## $var
## [1] 0.05833333
## 
## $df
## [1] 9
## 
## $a
## [1] 0.02759851
## 
## $b
## [1] 0.1944164

两个正态总体的置信区间

• 使用函数​​t.test()​​​进行$t$检验的一部分结果即为置信区间

均值差的置信区间

# 默认sigma未知，且不相等
interval_estimate2 <- function(x, y, 
   sigma=c(-1, -1), var.equal=FALSE, alpha=0.05) { 
   n1 <- length(x); 
   n2 <- length(y)
   xb <- mean(x); 
   yb <- mean(y)
   if (all(sigma >= 0)) 
{      tmp <- qnorm(1 - alpha / 2) * sqrt(sigma[1]^2 / n1 + sigma[2]^2 / n2)
      df <- n1 + n2
      }
   else {
      if (var.equal ==  TRUE) { 
        Sw <- ((n1 - 1)*var(x) + (n2 - 1)*var(y))/(n1 + n2 - 2)
        tmp <- sqrt(Sw*(1/n1 + 1/n2))*qt(1 - alpha/2,n1 + n2 - 2)
        df <- n1 + n2 - 2 
        }
      else {
        S1 <- var(x); 
        S2 <- var(y)
        nu <- (S1/n1 + S2/n2)^2 / (S1^2/n1^2/(n1 - 1) + S2^2/n2^2/(n2 - 1))
        tmp <- qt(1 - alpha/2, nu)*sqrt(S1/n1 + S2/n2)
        df <- nu
      }
     }
  list(mean=xb - yb, df=df, 
             a=xb - yb - tmp, b=xb - yb + tmp)
}

# 例题求解
# sigma未知时
set.seed(7)
x <- rnorm(100, 5.32, 2.18)
y <- rnorm(100, 5.76, 1.76)
interval_estimate2(x, y, sigma=c(2.18, 1.76))

## $mean
## [1] -0.3672189
## 
## $df
## [1] 200
## 
## $a
## [1] -0.9163587
## 
## $b
## [1] 0.1819209

set.seed(7)
x <- rnorm(12, 501.1, 2.4)
y <- rnorm(17, 499.7, 4.7)
interval_estimate2(x, y, var.equal=TRUE)

## $mean
## [1] 0.001928064
## 
## $df
## [1] 27
## 
## $a
## [1] -3.201143
## 
## $b
## [1] 3.204999

# 采用`t.test()`函数的方法
t.test(x, y, var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  x and y
## t = 0.0012351, df = 27, p-value = 0.999
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.201143  3.204999
## sample estimates:
## mean of x mean of y 
##  501.9227  501.9208

配对数据情形均值差的置信区间

配对数据作差，然后做单样本t检验，其中含有差的变化的区间估计

x <- c(11.3,15.0,15.0,13.5,12.8,10.0,11.0,12.0,13.0,12.3)
y <- c(14.0,13.8,14.0,13.5,13.5,12.0,14.7,11.4,13.8,12.0)
t.test(x-y)

## 
##  One Sample t-test
## 
## data:  x - y
## t = -1.3066, df = 9, p-value = 0.2237
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.8572881  0.4972881
## sample estimates:
## mean of x 
##     -0.68

方差比的区间估计

$\mu_1,\,\mu_2$已知

interval_var2 <- function(x, y, mu=c(Inf, Inf), alpha=0.05) { 
   n1 <- length(x); 
   n2 <- length(y) 
   # 均值已知
   if (all(mu < Inf)) {  
     Sx2<-1/n1*sum((x-mu[1])^2); 
     Sy2<-1/n2*sum((y-mu[2])^2)
     df1<-n1; 
     df2<-n2
   }
   # 均值未知
   else {      
     Sx2<-var(x); 
     Sy2<-var(y); 
     df1<-n1-1; 
     df2<-n2-1   
   }
   r <- Sx2/Sy2
   a <- r/qf(1-alpha/2,df1,df2)
   b <- r/qf(alpha/2,df1,df2)
   list(rate=r, df1=df1, df2=df2, a=a, b=b)
}

a <- c(79.98,80.04,80.02,80.04,80.03,80.03,80.04,79.97,80.05,80.03,80.02,80.00,80.02)
b <- c(80.02,79.94,79.98,79.97,79.97,80.03,79.95,79.97)
#均值已知μ1, μ2 =80
interval_var2(a, b, mu=c(80,80))

## $rate
## [1] 0.7326007
## 
## $df1
## [1] 13
## 
## $df2
## [1] 8
## 
## $a
## [1] 0.1760141
## 
## $b
## [1] 2.482042

#均值未知
interval_var2(a, b)

## $rate
## [1] 0.5837405
## 
## $df1
## [1] 12
## 
## $df2
## [1] 7
## 
## $a
## [1] 0.1251097
## 
## $b
## [1] 2.105269

非正态总体的区间估计

采用中心极限定理进行推导

首先进行数据标准化，当$n$充分大时，有
$\frac{\sum\limits_{i=1}^nX_i-u\mu}{\sqrt{n}\sigma}\sim N(0,\,1),$
参数$\mu$的区间估计($\sigma$已知)
$\left[\overline{X}-\frac{\sigma}{\sqrt n}Z_{\alpha/2},\,\overline{X}+\frac{\sigma}{\sqrt n}Z_{\alpha/2}\right]$

参数$\mu$的区间估计($\sigma$未知)

$\left[\overline{X}-\frac{S}{\sqrt n}Z_{\alpha/2},\,\overline{X}+\frac{S}{\sqrt n}Z_{\alpha/2}\right]$

编程得到

interval_estimate3<-function(x,sigma=-1,alpha=0.05) { 
   n<-length(x); 
   xb<-mean(x)
   if (sigma>=0)
      tmp<-sigma/sqrt(n)*qnorm(1-alpha/2)
   else
      tmp<-sd(x)/sqrt(n)*qnorm(1-alpha/2)
   list(mean=xb, a=xb-tmp, b=xb+tmp)
}

# 例题求解
x <- rexp(50,1/2.266)
interval_estimate3(x)

## $mean
## [1] 2.202523
## 
## $a
## [1] 1.654711
## 
## $b
## [1] 2.750334

单侧置信区间

单个总体均值的单侧置信区间

interval_estimate4<-function(x, sigma=-1, side=0, alpha=0.05){ 
   n<-length(x); xb<-mean(x)
   if (sigma>=0) { # σ已知
# side（标记），当标记<0时(左侧)，按置信上限公式求置信区间
      if (side<0) {         
        tmp<-sigma/sqrt(n)*qnorm(1-alpha)
        a <- -Inf; 
        b <- xb+tmp      
        }
      else if (side>0) {         
        tmp<-sigma/sqrt(n)*qnorm(1-alpha)
        a <- xb-tmp; 
        b <- Inf      
        }
      else {         
        tmp <- sigma/sqrt(n)*qnorm(1-alpha/2)
        a <- xb-tmp; b <- xb+tmp      
        } #默认side=0,求双侧置信区间
      df<-n   
      }
   else {      
     if (side<0) {         
       tmp <- sd(x)/sqrt(n)*qt(1-alpha,n-1)
         a <- -Inf; 
         b <- xb+tmp      
         }
      else if (side>0) {        
        tmp <- sd(x)/sqrt(n)*qt(1-alpha,n-1)
        a <- xb-tmp; b <- Inf      
        }
      else {         
        tmp <- sd(x)/sqrt(n)*qt(1-alpha/2,n-1)
         a <- xb-tmp; b <- xb+tmp      
         } #求双侧置信区间
      df<-n-1   
      }
   list(mean=xb, df=df, a=a, b=b)
}

# 例题求解
x <- c(1050,1100,1120,1250,1280)
interval_estimate4(x, side=1)

## $mean
## [1] 1160
## 
## $df
## [1] 4
## 
## $a
## [1] 1064.9
## 
## $b
## [1] Inf

单个总体方差的单侧置信区间

interval_var3<-function(x,mu=Inf,side=0,alpha=0.05) { 
   n<-length(x)
   if (mu<Inf) {      
     S2<-sum((x-mu)^2)/n; df<-n   
     }
   else {      
     S2<-var(x); df<-n-1   
     }
   if (side<0) {      
     a <- 0
     b <- df*S2/qchisq(alpha,df)   
      }
   else if (side>0) {      
     a <- df*S2/qchisq(1-alpha,df)
     b <- Inf   
   }
   else {      
     a<-df*S2/qchisq(1-alpha/2,df) 
     b<-df*S2/qchisq(alpha/2,df)   
     }
list(var=S2, df=df, a=a, b=b)
}

# 例题求解
x <- c(10.1,10,9.8,10.5,9.7,10.1,9.9,10.2,10.3,9.9)
interval_var3(x, side=-1)

## $var
## [1] 0.05833333
## 
## $df
## [1] 9
## 
## $a
## [1] 0
## 
## $b
## [1] 0.1578894

两个总体均值差的单侧置信区间

interval_estimate5<-function(x, y,sigma=c(-1,-1), var.equal=FALSE, side=0, alpha=0.05) {
  n1<-length(x); n2<-length(y)
  xb<-mean(x); yb<-mean(y); zb<-xb-yb
  if (all(sigma>=0)){
    if (side<0){
    tmp<-qnorm(1-alpha)*sqrt(sigma[1]^2/n1+sigma[2]^2/n2)
    a <- -Inf; b <- zb+tmp
    }
  else if (side>0){
    tmp<-qnorm(1-alpha)*sqrt(sigma[1]^2/n1+sigma[2]^2/n2)
    a <- zb-tmp; b <- Inf
  }
  else{
    tmp<-qnorm(1-alpha/2)*sqrt(sigma[1]^2/n1+sigma[2]^2/n2)
    a <- zb-tmp; b <- zb+tmp
  }
  df<-n1+n2
  }
  else{
    if (var.equal == TRUE){
      Sw<-((n1-1)*var(x)+(n2-1)*var(y))/(n1+n2-2)
      if (side<0){
        tmp<-sqrt(Sw*(1/n1+1/n2))*qt(1-alpha,n1+n2-2)
        a <- -Inf; b <- zb+tmp
      }
    else if (side>0){
      tmp<-sqrt(Sw*(1/n1+1/n2))*qt(1-alpha,n1+n2-2)
      a <- zb-tmp; b <- Inf
    }
    else{
      tmp<-sqrt(Sw*(1/n1+1/n2))*qt(1-alpha/2,n1+n2-2)
      a <- zb-tmp; b <- zb+tmp
    }
    df<-n1+n2-2
  }
  else{
    S1<-var(x); S2<-var(y)
    nu<-(S1/n1+S2/n2)^2/(S1^2/n1^2/(n1-1)+S2^2/n2^2/(n2-1))
      if (side<0){
        tmp<-qt(1-alpha, nu)*sqrt(S1/n1+S2/n2)
        a <- -Inf; b <- zb+tmp
      }
      else if (side>0){
        tmp<-qt(1-alpha, nu)*sqrt(S1/n1+S2/n2)
        a <- zb-tmp; b <- Inf
      }
      else{
        tmp<-qt(1-alpha/2, nu)*sqrt(S1/n1+S2/n2)
        a <- zb-tmp; b <- zb+tmp
      }
    df<-nu
    }
  }
list(mean=zb, df=df, a=a, b=b)
}

两个总体方差的置信区间

interval_var4<-function(x,y,mu=c(Inf, Inf), side=0, alpha=0.05){
  n1<-length(x); n2<-length(y)
  if (all(mu<Inf)) {
    Sx2<-1/n1*sum((x-mu[1])^2); df1<-n1
    Sy2<-1/n2*sum((y-mu[2])^2); df2<-n2
  }
  else{
    Sx2<-var(x); Sy2<-var(y); df1<-n1-1; df2<-n2-1
  }
  r<-Sx2/Sy2
  if (side<0) {
    a <- 0
    b <- r/qf(alpha,df1,df2)
  }
  else if (side>0) {
    a <- r/qf(1-alpha,df1,df2)
    b <- Inf
  }
  else{
    a<-r/qf(1-alpha/2,df1,df2)
    b<-r/qf(alpha/2,df1,df2)
  }
  list(rate=r, df1=df1, df2=df2, a=a, b=b)
}