R语言怎么算估计的绝对误差 r语言计算矩估计

第五章 参数估计

  学习《R语言与统计分析》（汤银才 主编）一书，整理第五章参数估计部分如下：

5.1 矩法估计和极大似然估计

5.1.1 矩法估计

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

例：

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

5.1.2 极大似然估计

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

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

例：

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

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

5.2.1 均值$\mu$的区间估计

1. 方差$\sigma^2$已知时$\mu$的置信区间

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

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)
}

验证一个例子

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()

例：

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

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

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

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)
}

给出一个验证

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

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

5.3.1 均值差$\mu_1-\mu_2$的置信区间

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

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

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)
}

看一个例子：

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)

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

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

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)

5.3.2 均值差$\sigma^2_1/\sigma^2_2$的置信区间

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

例子：

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)

5.4 单总体比率$p$的区间估计

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

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

一个例子 ：矫正

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

不矫正：

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

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

binom.test(38, 200)

5.5 两总体比率差$p_i-p_2$的区间估计

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

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

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

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)
}

给一个例子

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

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

1. 总体方差$\sigma^2$已知

给出自定义函数size.norm1：

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，应抽取多少户进行调查

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

2. 总体方差$\sigma^2$未知

给出函数size.norm2()

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=10$，试问$\alpha=0.01$下要抽多少样本

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

5.6.2 估计比例$p$时样本容量的确定

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

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

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

size.bin(0.03, 0.9, 0.95)