写在前面

这次总结一下数理统计中的参数估计，即**点估计(矩估计、极大似然估计)和区间估计(置信区间)**部分的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;θ),θ∈Θf(x;\,\theta),\,\theta\in\Thetaf(x;θ),θ∈Θ是未知参数，X1,X2,⋯,XnX_1,\,X_2,\,\cdots,\,X_nX1​,X2​,⋯,Xn​为来自总体XXX的样本，称
  L(θ;x)=L(θ;x1,x2,⋯,xn)=∏i=1nf(xi;θ)L(\theta;\,x)=L(\theta;x_1,\,x_2,\,\cdots,\,x_n)=\prod\limits_{i=1}^nf(x_i;\,\theta) L(θ;x)=L(θ;x1​,x2​,⋯,xn​)=i=1∏n​f(xi​;θ)
  为θ\thetaθ的极大似然函数(likelihood function)。

• 定义：设总体X的概率密度函数或分布律为f(x;θ),θ∈Θf(x;\,\theta),\,\theta\in\Thetaf(x;θ),θ∈Θ是未知参数，X1,X2,⋯,XnX_1,\,X_2,\,\cdots,\,X_nX1​,X2​,⋯,Xn​为来自总体XXX的样本，L(θ;x)L(\theta;\,x)L(θ;x)为θ\thetaθ的似然函数, 若θ^=θ^(X)=θ^(X1,X2,⋯,Xn)\hat{\theta}=\hat{\theta}(X)=\hat{\theta}(X_1,\,X_2,\,\cdots,\,X_n)θ^=θ^(X)=θ^(X1​,X2​,⋯,Xn​)是一个统计量，且满足：
  L(θ^(X);X)=sup⁡θ∈ΘL(θ;X)L(\hat{\theta}(X);\,X)=\sup\limits_{\theta\in\Theta}L(\theta;\,X) L(θ^(X);X)=θ∈Θsup​L(θ;X)
  则称θ^\hat{\theta}θ^为θ\thetaθ的最大似然估计。

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

• 均匀分布

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

• 指数分布

  服从指数分布的最大似然估计函数为
  L(λ;x)=λne−λ∑i=1nxiL(\lambda;\,x) =\lambda^n\mathrm{e}^{-\lambda\sum\limits_{i=1}^nx_i} L(λ;x)=λne−λi=1∑n​xi​
  取对数并求导得到
  ∂ln⁡L(λ;x)∂λ=(nln⁡λ−λ∑i=1nxi)λ=nλ−∑i=1nxi=0\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 ∂λ∂lnL(λ;x)​=(nlnλ−λi=1∑n​xi​)λ​=λn​−i=1∑n​xi​=0
  即λ=n∑i=1nxi\lambda=\dfrac{n}{\sum\limits_{i=1}^nx_i}λ=i=1∑n​xi​n​.

• 正态分布

  正态分布的似然函数为

L(μ,σ2;x)=∏i=1nf(xi;μ,σ2)=(2πσ2)−n2exp⁡[−12σ2∑i=1n(xi−μ)2],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], L(μ,σ2;x)=i=1∏n​f(xi​;μ,σ2)=(2πσ2)−2n​exp[−2σ21​i=1∑n​(xi​−μ)2],

对数似然函数为
ln⁡L(μ,σ2;x)=−n2ln⁡(2πσ2)−12σ2∑i=1n(xi−μ)2,\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, lnL(μ,σ2;x)=−2n​ln(2πσ2)−2σ21​i=1∑n​(xi​−μ)2,
令
{∂ln⁡L(μ,σ2;x)∂μ=1σ2∑i=1n(xi−μ)=0,∂ln⁡L(μ,σ2;x)∂σ2=−n2σ2+12σ4∑i=1n(xi−μ)2=0,\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} ⎩⎪⎪⎨⎪⎪⎧​∂μ∂lnL(μ,σ2;x)​=σ21​i=1∑n​(xi​−μ)=0,∂σ2∂lnL(μ,σ2;x)​=−2σ2n​+2σ41​i=1∑n​(xi​−μ)2=0,​
解此似然方程组得到：
μ=1n∑i=1nxi=x‾,σ2=1n∑i=1n(xi−x‾)2,\mu=\dfrac1n\sum\limits_{i=1}^nx_i=\overline{x},\quad \sigma^2=\dfrac1n\sum_{i=1}^n(x_i-\overline{x})^2, μ=n1​i=1∑n​xi​=x,σ2=n1​i=1∑n​(xi​−x)2,
进一步验证，对于对数似然函数ln⁡L(μ,σ2;x)\ln L(\mu,\,\sigma^2;\,x)lnL(μ,σ2;x)的二阶Hesse矩阵
[−nσ2−1σ4∑i=1n(xi−μ)−1σ4∑i=1n(xi−μ)n2σ4−1σ6∑i=1n(xi−μ)2]=[−nσ200−n2σ4]\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} ⎣⎢⎡​−σ2n​−σ41​i=1∑n​(xi​−μ)​−σ41​i=1∑n​(xi​−μ)2σ4n​−σ61​i=1∑n​(xi​−μ)2​⎦⎥⎤​=⎣⎡​−σ2n​0​0−2σ4n​​⎦⎤​
为负定矩阵，所以(x‾,1n∑i=1n(xi−x‾)2)\left(\overline{x},\,\dfrac1n\sum\limits_{i=1}^n(x_i-\overline{x})^2\right)(x,n1​i=1∑n​(xi​−x)2)是L(μ,σ2;x)L(\mu,\,\sigma^2;\,x)L(μ,σ2;x)的极大值点。故(μ,σ2)(\mu,\,\sigma^2)(μ,σ2)的最大似然估计为
μ^=X‾,σ^2=1n∑i=1n(Xi−X‾)2.\hat{\mu}=\overline{X},\quad\hat{\sigma}^2=\dfrac1n\sum\limits_{i=1}^n(X_i-\overline{X})^2. μ^​=X,σ^2=n1​i=1∑n​(Xi​−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()函数，求解有nnn个方程、nnn个未知量的非线性方程组。

# 定义待求方程
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=∫0θx1θdx=θ2,EX=\int_0^\theta x\frac1\theta\mathrm{d}x=\frac\theta2, EX=∫0θ​xθ1​dx=2θ​,
    所以其矩估计为θ^=2X‾=2n∑i=1nXi\hat{\theta}=2\overline{X}=\dfrac2n\sum\limits_{i=1}^nX_iθ^=2X=n2​i=1∑n​Xi​.

  • 情形二（一般情况）

  EX=∫abx1b−adx=b+a2,DX=∫abx21b−adx−(b+a2)2=(b−a)212,\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} EXDX​=∫ab​xb−a1​dx=2b+a​,=∫ab​x2b−a1​dx−(2b+a​)2=12(b−a)2​,​
  令
  {b+a2=X‾(b−a)212=1n∑i=1nXi2\begin{cases} \dfrac{b+a}2=\overline{X}\\ \dfrac{(b-a)^2}{12}=\dfrac1n\sum\limits_{i=1}^nX_i^2 \end{cases} ⎩⎪⎨⎪⎧​2b+a​=X12(b−a)2​=n1​i=1∑n​Xi2​​
  解得a^=X‾−3n∑i=1nXi2,b^=X‾+3n∑i=1nXi2.\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}.a^=X−n3​i=1∑n​Xi2​​, b^=X+n3​i=1∑n​Xi2​​.

• 指数分布
  EX=∫0+∞λxe−λxdx=1λ,EX=\int_0^{+\infty}\lambda x\mathrm{e}^{-\lambda x}\mathrm{d}x=\frac1\lambda, EX=∫0+∞​λxe−λxdx=λ1​,
  因此其矩估计为λ^=n∑i=1nXi\hat{\lambda}=\dfrac{n}{\sum\limits_{i=1}^{n}X_i}λ^=i=1∑n​Xi​n​.

• 正态分布

  算总体XXX的一阶、二阶原点矩

  M1=EX=μ,M2=EX2=σ2+μ2M_1 =EX=\mu,\quad M_2 =EX^2=\sigma^2+\mu^2 M1​=EX=μ,M2​=EX2=σ2+μ2
  以及样本的一阶、二阶原点矩

  A1=X‾=1n∑i=1nXi,A2=1n∑i=1nXi2.A_1=\overline{X}=\frac1n\sum_{i=1}^nX_i,\quad A_2=\frac1n\sum_{i=1}^nX_i^2. A1​=X=n1​i=1∑n​Xi​,A2​=n1​i=1∑n​Xi2​.

  所以得到方程组

  {μ=X‾σ2+μ2=1n∑i=1nXi2\begin{cases} \mu=\overline{X}\\ \sigma^2+\mu^2=\dfrac1n\sum\limits_{i=1}^nX_i^2 \end{cases} ⎩⎨⎧​μ=Xσ2+μ2=n1​i=1∑n​Xi2​​

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

  μ^=X‾,σ^2=1n∑i=1nXi2−X‾2=1n∑i=1n(Xi−X‾)2=n−1nS2.\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} μ^​σ^2​=X,=n1​i=1∑n​Xi2​−X2=n1​i=1∑n​(Xi​−X)2=nn−1​S2.​

# 定义待求方程组
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条件分支语句。

一个正态总体的置信区间

σ2\sigma^2σ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

σ2\sigma^2σ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

方差σ2\sigma^2σ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()进行ttt检验的一部分结果即为置信区间

均值差的置信区间

# 默认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

方差比的区间估计

μ1,μ2\mu_1,\,\mu_2μ1​,μ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/Sy2a <- 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

非正态总体的区间估计

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

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

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

[X‾−SnZα/2,X‾+SnZα/2]\left[\overline{X}-\frac{S}{\sqrt n}Z_{\alpha/2},\,\overline{X}+\frac{S}{\sqrt n}Z_{\alpha/2}\right] [X−n​S​Zα/2​,X+n​S​Zα/2​]

编程得到

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)elsetmp<-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 <- 0b <- 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-ybif (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<-n1Sy2<-1/n2*sum((y-mu[2])^2); df2<-n2}else{Sx2<-var(x); Sy2<-var(y); df1<-n1-1; df2<-n2-1}r<-Sx2/Sy2if (side<0) {a <- 0b <- 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)
}

R语言学习笔记（四）参数估计相关推荐

1. R语言学习笔记4_参数估计

   目录 四.参数估计 4.1 矩估计和极大似然估计法 4.1.1 矩估计 4.1.2 极大似然估计 单参数 optimize( ) 多参数 optim( ) .nlm( ) 4.2 单正态总体参数的区间 ...

2. R语言学习笔记——高级篇：第十四章-主成分分析和因子分析

   R语言 R语言学习笔记--高级篇:第十四章-主成分分析和因子分析 文章目录 R语言 前言 一.R中的主成分和因子分析 二.主成分分析 2.1.判断主成分的个数 2.2.提取主成分 2.3.主成分旋转 ...

3. R语言学习笔记（1~3）

   R语言学习笔记(1~3) 一.R语言介绍 x <- rnorm(5) 创建了一个名为x的向量对象,它包含5个来自标准正态分布的随机偏差. 1.1 注释 由符号#开头. #函数c()以向量的形式输 ...

4. R语言学习笔记——入门篇：第一章-R语言介绍

   R语言 R语言学习笔记--入门篇:第一章-R语言介绍 文章目录 R语言 一.R语言简介 1.1.R语言的应用方向 1.2.R语言的特点 二.R软件的安装 2.1.Windows/Mac 2.2.Lin ...

5. R语言学习笔记——入门篇：第三章-图形初阶

   R语言 R语言学习笔记--入门篇:第三章-图形初阶 文章目录 R语言 一.使用图形 1.1.基础绘图函数:plot( ) 1.2.图形控制函数:dev( ) 补充--直方图函数:hist( ) 补充- ...

6. r语言c函数怎么用,R语言学习笔记——C#中如何使用R语言setwd()函数

   在R语言编译器中,设置当前工作文件夹可以用setwd()函数. > setwd("e://桌面//") > setwd("e:\桌面\") > ...

7. R语言学习笔记 07 Probit、Logistic回归

   R语言学习笔记 文章目录 R语言学习笔记 probit回归 factor()和as.factor() relevel() 案例11.4复刻 glm函数 整理变量 回归:Logistic和Probit- ...

8. R语言学习笔记 06 岭回归、lasso回归

   R语言学习笔记 文章目录 R语言学习笔记 比较lm.ridge和glmnet函数 画岭迹图 图6-4 <统计学习导论 基于R语言的应用>P182 图6-6<统计学习导论 基于R语言的 ...

9. R语言学习笔记（八）--读写文件与网络爬虫

   R语言学习笔记(八) 1 工作路径 2 保存R对象 3 Scan函数 3-1 从控制台读取数据 3-2 从txt文件读取数据 3-3 从url读取数据 4 按行读写文本文件 5 读取文本文件(txt. ...

10. R语言学习笔记（三）多元数据的数据特征、相关分析与图形表示

    文章目录 写在前面 独立性检验 χ2\chi^2χ2独立性检验 Fisher独立性检验 Cochran-Mantel-Haenszel χ2\chi^2χ2独立性检验 相关性分析 相关性检验 相关性检 ...

最新文章

1. 你竟然还不懂变分自编码机？这个16岁的OpenAI天才实习生讲得可透彻了
2. 我们团队设计的一个基于微服务的高并发服务器架构
3. CSS 如何让Table的里面TD全有边框 而Table的右左边框没有
4. 前后端数据加密传输 RSA非对称加密
5. 【发现问题】Java中PrintStream和PrintWriter的区别
6. 分段线性变换与直方图修正
7. liunx 命令手册 （chm）
8. ECCV18 | 无监督难分样本挖掘改进目标检测
9. 树莓派智能家居-语音聊天机器人实现
10. [代码发布]中文文字转换组件 1.0,支持VB/ASP编程
11. 高一信息技术 计算机配件的真伪辨别,高一信息技术组PPT.ppt
12. c语言编写ocr软件,开源OCR引擎Tesseract
13. QT界面程序异常结束问题分析 ，弹出 SogouInput\Components\
14. java把字符串转为日期_Java程序将字符串转换为日期
15. 从appfuse开始学习Spring和Hibernate - （1）构建项目
16. 渗透之——触发Easy File Sharing Web Server 7.2 HEAD缓冲区溢出的Python脚本
17. 源码解析kafka删除topic
18. Unity SteamVR报错问题却影响运行的记录（Log path could not be located (112)“）
19. 微信小程序 内容换行
20. 【Pyecharts50例】自定义饼图标签/显示百分比

热门文章

1. EXCEL下载功能在XP系统中运行是好好的,到windows2003系统上,就报错
2. C语言指针作为参数的传递问题
3. 广度优先搜索——Corn Maze S(洛谷 P1825)
4. 贪心算法—区间调度 电影节(POJ 4151)
5. java文件流下载excel_React获取Java后台文件流下载Excel文件
6. unity 自动将文件上传_unity如何存储文件夹
7. sqlplus登录缓慢的问题分析过程及解决小记
8. PostgreSQL实际场景的十大缺陷你知道吗？
9. SQL执行效率提升几万倍的操作详解！
10. 推荐几个阿里、美团、腾讯大佬的公众号，一起学习！