更新于 

参数估计

矩法估计和极大似然估计

矩法估计

无固定的函数，矩估计需自行构造统计量，并计算相关数据

例：

1
2
3
4
X<-sample(c(1,0), 20, replace = T, prob = c(0.6, 0.4))
theta <- mean(X)
t <- theta/(1-theta)  # 构造的t统计量
t

在这里插入图片描述

极大似然估计

需自行求解出(对数)极大似然函数，可以使用optimize()函数求解极大值和极大值点。多参数时可使用optim()或nlm()

optimize(f = , interval = ,lower = min(interval), upper = max(interval), maximum = T, tol = .Machine$dobule.eps^{0.25}, \cdots

例：

1
2
f <- function(p)(p^517)*(1-p)^483
optimize(f, c(0,1), maximum = T)

在这里插入图片描述

单正态总体参数的区间估计

均值μ\mu的区间估计

方差σ2\sigma^2已知时μ\mu的置信区间

这种情况下R语言中没有相关内置函数，以下给出方差σ2\sigma^2已知时μ\mu的置信区间估计的R语言函数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
z.test <- function(x, n, sigma, alpha=0.05, u0=0, alternative="two.side"){
  mean <- mean(x)
  z <- (mean-u0)/(sigma/sqrt(n))
  p <- pnorm(z, lower.tail = F)
  mean <- mean
  z <- z
  p.value <- p
  if(alternative =="two.side"){
    p.value <- 2*pnorm(abs(z), lower.tail = F)
  }else{
    p.value <- pnorm(z)
  }
  conf.int <- c(
    mean-sigma*qnorm(1-alpha/2, mean=0, sd=1, 
                     lower.tail = T)/sqrt(n),
    mean+sigma*qnorm(1-alpha/2, mean=0, sd=1, 
                     lower.tail = T)/sqrt(n)
  )
  out <- list(mean=mean, z=z, p.value=p.value, conf.int=conf.int)
}

验证一个例子

1
2
3
x <- c(175, 176, 173, 175, 174, 173, 173, 176, 173, 179)
result <-z.test(x, 10, 1.5, 0.05)
result

在这里插入图片描述

方差σ2\sigma^2未知时μ\mu的置信区间

此时R语言中有内置函数t.test()

例：

1
2
3
x <- c(175, 176, 173, 175, 174, 173, 173, 176, 173, 179)
result <-t.test(x)
result

在这里插入图片描述

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

讨论μ\mu未知，此种情况R语言中依然没有内置函数求解方差σ2\sigma^2的区间估计，给出函数chisq.var.test如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
chisq.var.test <- function(x, var, alpha, alternative = "two.side"){
  n <- length(x)
  v <- var(x)
  chi <- (n-1)*v/var
  p.value <- pchisq(chi, n-1)
  
  if(alternative=="less"){
    p.value <- pchisq(chi, n-1, lower.tail = F)
  }else if(alternative=="two.side"){
    p.value <- 2*min(pchisq(chi, n-1), 
                     pchisq(chi, n-1, lower.tail = F))
  }
  
  conf.int <- c(
    (n-1)*v/qchisq(alpha/2, df=n-1, lower.tail = F),
    (n-1)*v/qchisq(alpha/2, df=n-1, lower.tail = T)
  )
  
  out <- list(var=v, chi2=chi, p.value=p.value, conf.int=conf.int)
}

给出一个验证

1
2
3
x <- c(175, 176, 173, 175, 174, 173, 173, 176, 173, 179)
result <- chisq.var.test(x , var=1.5, alpha=0.05)
result

在这里插入图片描述

正态总体参数的区间估计

均值差μ1−μ2\mu_1-\mu_2的置信区间

两方差都已知时两均值差的置信区间

依然没有内置函数，需给出函数two.sample.ci()，求解双侧置信区间

1
2
3
4
5
6
7
two.sample.ci <- function(x, y, conf.int = 0.95, sigma1, sagma2){
  alpha = 1-conf.int
  m <- length(x); n<-length(y)
  xbar = mean(x)-mean(y)
  zstar = qnorm(1-alpha/2)*sqrt((sigma1/m+sagma2/n))
  xbar + c(-zstar,+zstar)
}

看一个例子：

1
2
3
4
x <- c(628, 583, 510, 554, 612, 523, 530, 615)
y <- c(534, 433, 398, 470, 567, 480, 498, 560, 503, 426)
sigma1 <- 2140;sigma2<-3250
two.sample.ci(x, y, conf.int=0.95, sigma1, sigma2)

在这里插入图片描述

两方差都未知但相等时两均值差的置信区间

此时可以直接利用t.test()函数

1
2
3
x <- c(628, 583, 510, 554, 612, 523, 530, 615)
y <- c(534, 433, 398, 470, 567, 480, 498, 560, 503, 426)
t.test(x,y,var.equal = T)

在这里插入图片描述

均值差σ12/σ22\sigma^2_1/\sigma^2_2的置信区间

此时,在R中var.test()可以直接用于求两正态总体方差比的置信区间

例子：

1
2
3
x <- c(20.5, 19.8, 19.7, 20.4, 20.1, 20.0, 19.0, 19.9)
y <- c(20.7, 19.8, 19.5, 20.8, 20.4, 19.6, 20.2)
var.test(x,y)

在这里插入图片描述

单总体比率pp的区间估计

假设估计在总体中具有某种特性的个体站总体的比例，设为pp。例如整个学校中女生占全校人数的比例，产品的不合格率，电视的收视率，政策的支持率等。这里采用一种求pp的近似区间估计的方法，
称在样本中具有某种特征的个体站样本总数的比例为样本比例。设xx为容量为n的样本中具有某种特征的个体数量，则样本比例为x/nx/n。当总体中的样品数足够多时，xx近似服从二项分布b(n,p)b(n,p)(实际上它是超几何分布)，这时总体比例可用样本比例来估计，及p^=xn\hat p=\frac{x}{n}，且为极大似然估计。

可直接使用prop.test()对pp进行估计与检验，这时用正态分布来近似

一个例子 ：矫正

1
prop.test(38, 200, correct = T)

在这里插入图片描述

不矫正：

1
prop.test(38, 200, correct = F)

在这里插入图片描述

用二项分布来近似binom.test()

1
binom.test(38, 200)

在这里插入图片描述

两总体比率差pi−p2p_i-p_2的区间估计

这里仍可使用函数prop.test()，这里给出的是经过连续性修正之后的结果

1
2
3
like <- c(478, 246)
people <- c(1000, 750)
prop.test(like, people)

在这里插入图片描述

可自行给出函数ratio.ci() 计算没有修正的两比例之间的区间估计

1
2
3
4
5
6
7
8
ratio.ci <- function(x,y,n1,n2,conf.level){
  xbar1 <- x/n1; xbar2 <- y/n2
  xbar <- xbar1-xbar2
  alpha <- 1-conf.level
  zstar <- qnorm(1-alpha/2)*
    (xbar1*(1-xbar1)/n1+xbar2*(1-xbar2)/n2)^(1/2)
  xbar+c(-zstar,+zstar)
}

给一个例子

1
ratio.ci(478, 246, 1000, 750, conf.level = 0.95)

在这里插入图片描述

估计正态总体均值时样本容量的确定

总体方差σ2\sigma^2已知

给出自定义函数size.norm1：

1
2
3
4
size.norm1 <- function(d, var, conf.level=0.95){
  alpha <- 1-conf.level
  ((qnorm(1-alpha/2)*var^(1/2))/d)^2
}

给一个例子：某地区1000户，拟抽取一个简单的样本调查一个月的评价开始，要求置信水平95%，允许最大误差为2，根据经验，家庭间开支的方差为500，应抽取多少户进行调查

1
size.norm1(2, 500, conf.level = 0.95)

在这里插入图片描述

总体方差σ2\sigma^2未知

给出函数size.norm2()

1
2
3
4
5
6
7
8
9
10
11
size.norm2 <- function(s, alpha, d, m){
  t0 <- qt(alpha/2, m, lower.tail = F)
  n0 <- (t0*s/d)^2
  t1 <- qt(alpha/2, n0, lower.tail = F)
  n1 <- (t1*s/d)^2
  while(abs(n1-n0)>0.5){
    n0 <- (qt(alpha/2, n1, lower.tail = F)*s/d)^2
    n1 <- (qt(alpha/2, n0, lower.tail = F)*s/d)^2
  }
  n1
}

给一个例子：

某公司生产了一批新产品，产品总体服从正态分布，现要估计这批产品的平均重量，最大允许误差2，样本标准差s=10s=10，试问α=0.01\alpha=0.01下要抽多少样本

1
size.norm2(10, 0.01, 2, 100)

在这里插入图片描述

估计比例pp时样本容量的确定

同样给出自定义函数size.bin <- f()

1
2
3
4
size.bin <- function(d, p, conf.level = 0.95){
  alpha = 1 - conf.level
  ((qnorm(1-alpha/2))/d)^2*p*(1-p)
}

例 ：某市一所重点大学历届毕业生就业率为90%，估计应用毕业生就业率，要求估计误差不超过3%，试问α=0.05\alpha=0.05下要抽取应届毕业生多少人？

1
size.bin(0.03, 0.9, 0.95)

在这里插入图片描述