抱歉,您的浏览器无法访问本站

本页面需要浏览器支持(启用)JavaScript


了解详情 >

  这是曾经R语言软件课程的简单笔记的参数估计部份,整理至此。

矩法估计和极大似然估计

矩法估计

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

例:

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)

在这里插入图片描述

两总体比率差pip2p_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)

在这里插入图片描述




萌ICP备20201324号 | 渝ICP备2021006569号 | 公安备50011702500564号

又拍云logo     腾讯云logo