文章目录

写在前面点估计极大似然估计可求出解析解不易或无法求出解析解 矩估计 区间估计一个正态总体的置信区间${\sigma }^2$已知时，$\mu$的区间估计${\sigma }^2$未知时，$\mu$的区间估计方差${\sigma }^2$的区间估计 两个正态总体的置信区间均值差的置信区间配对数据情形均值差的置信区间 方差比的区间估计 非正态总体的区间估计单侧置信区间单个总体均值的单侧置信区间单个总体方差的单侧置信区间两个总体均值差的单侧置信区间两个总体方差的置信区间

写在前面

这次总结一下数理统计中的参数估计，即**点估计(矩估计、极大似然估计)和区间估计(置信区间)**部分的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 _{i=1}^{n}f(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)=\underset{\theta \in \Theta }{\sup }L(\theta ;X)$$ 则称$\hat{\theta }$为$\theta$的最大似然估计。

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

均匀分布

显然得到$\hat{a}=X_{(1)},\hat{b}=X_{(n)}$.

指数分布

服从指数分布的最大似然估计函数为 $$L(\lambda ;x)={\lambda }^n{\mathrm{e}}^{-\lambda \sum _{i=1}^n{x}_i}$$ 取对数并求导得到 $$\frac{\partial \ln L(\lambda ;x)}{\partial \lambda }={\left(n\ln \lambda -\lambda \sum _{i=1}^n{x}_i\right)}_{\lambda }=\frac{n}{\lambda }-\sum _{i=1}^n{x}_i=0$$ 即$\lambda =\frac{n}{\sum _{i=1}^n{x}_i}$.

正态分布

正态分布的似然函数为

$$L(\mu ,{\sigma }^2;x)=\prod _{i=1}^{n}f(x_i;\mu ,{\sigma }^2)=(2\pi {\sigma }^2{)}^{-\frac{n}{2}}\exp \left[-\frac{1}{2{\sigma }^2}\sum _{i=1}^n(x_i-\mu {)}^2\right],$$

对数似然函数为 $$\ln L(\mu ,{\sigma }^2;x)=-\frac{n}{2}\ln (2\pi {\sigma }^2)-\frac{1}{2{\sigma }^2}\sum _{i=1}^n(x_i-\mu {)}^2,$$ 令 $$\{\begin{array}{l}\frac{\partial \ln L(\mu ,{\sigma }^2;x)}{\partial \mu }=\frac{1}{{\sigma }^2}\sum _{i=1}^n(x_i-\mu )=0,\\ \frac{\partial \ln L(\mu ,{\sigma }^2;x)}{\partial {\sigma }^2}=-\frac{n}{2{\sigma }^2}+\frac{1}{2{\sigma }^4}\sum _{i=1}^n(x_i-\mu {)}^2=0,\end{array}$$ 解此似然方程组得到： $$\mu =\frac{1}{n}\sum _{i=1}^n{x}_i=\overline{x},{\sigma }^2=\frac{1}{n}\sum _{i=1}^n(x_i-\overline{x}{)}^2,$$ 进一步验证，对于对数似然函数$\ln L(\mu ,{\sigma }^2;x)$的二阶Hesse矩阵 $$\left[\begin{array}{cc}-\frac{n}{{\sigma }^2}& -\frac{1}{{\sigma }^4}\sum _{i=1}^n(x_i-\mu )\\ -\frac{1}{{\sigma }^4}\sum _{i=1}^n(x_i-\mu )& \frac{n}{2{\sigma }^4}-\frac{1}{{\sigma }^6}\sum _{i=1}^n(x_i-\mu {)}^2\end{array}\right]=\left[\begin{array}{cc}-\frac{n}{{\sigma }^2}& 0\\ 0& -\frac{n}{2{\sigma }^4}\end{array}\right]$$ 为负定矩阵，所以$\left(\overline{x},\frac{1}{n}\sum _{i=1}^n(x_i-\overline{x}{)}^2\right)$是$L(\mu ,{\sigma }^2;x)$的极大值点。故$(\mu ,{\sigma }^2)$的最大似然估计为 $$\hat{\mu }=\overline{X},{\hat{\sigma }}^2=\frac{1}{n}\sum _{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\frac{1}{\theta }\mathrm{d}x=\frac{\theta }{2},$$ 所以其矩估计为$\hat{\theta }=2\overline{X}=\frac{2}{n}\sum _{i=1}^n{X}_i$.

情形二（一般情况）

$$\begin{array}{rl}EX& ={\int }_a^{b}x\frac{1}{b-a}\mathrm{d}x=\frac{b+a}{2},\\ DX& ={\int }_a^b{x}^2\frac{1}{b-a}\mathrm{d}x-{\left(\frac{b+a}{2}\right)}^2=\frac{(b-a{)}^2}{12},\end{array}$$ 令 $$\{\begin{array}{l}\frac{b+a}{2}=\overline{X}\\ \frac{(b-a{)}^2}{12}=\frac{1}{n}\sum _{i=1}^n{X}_i^2\end{array}$$ 解得$\hat{a}=\overline{X}-\sqrt{\frac{3}{n}\sum _{i=1}^n{X}_i^2},\hat{b}=\overline{X}+\sqrt{\frac{3}{n}\sum _{i=1}^n{X}_i^2}.$

指数分布 $$EX={\int }_0^{+\infty }\lambda x{\mathrm{e}}^{-\lambda x}\mathrm{d}x=\frac{1}{\lambda },$$ 因此其矩估计为$\hat{\lambda }=\frac{n}{\sum _{i=1}^n{X}_i}$.

正态分布

算总体$X$的一阶、二阶原点矩

$$M_1=EX=\mu ,M_2=E{X}^2={\sigma }^2+{\mu }^2$$ 以及样本的一阶、二阶原点矩

$$A_1=\overline{X}=\frac{1}{n}\sum _{i=1}^n{X}_i,A_2=\frac{1}{n}\sum _{i=1}^n{X}_i^2.$$

所以得到方程组

$$\{\begin{array}{l}\mu =\overline{X}\\ {\sigma }^2+{\mu }^2=\frac{1}{n}\sum _{i=1}^n{X}_i^2\end{array}$$

解上述方程，得均值$\mu$和方差${\sigma }^2$的矩估计

$$\begin{array}{rl}\hat{\mu }& =\overline{X},\\ {\hat{\sigma }}^2& =\frac{1}{n}\sum _{i=1}^n{X}_i^2-{\overline{X}}^2=\frac{1}{n}\sum _{i=1}^n(X_i-\overline{X}{)}^2=\frac{n-1}{n}{S}^2.\end{array}$$

# 定义待求方程组 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 _{i=1}^n{X}_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) }