2

我试图找到给定 beta = 1 的 beta 分布的 alpha 的 MLE 估计值。我尝试使用来自estimatetools 包的 maxlogL 但 g

x <- rbeta(n = 1000, shape1 = 0.7, shape2 = 1)

alpha_hat <- maxlogL(x = x, dist = "dbeta", fixed = list(shape2 = 1), lower = (0), upper = (1), link = list(over = "shape1", fun = "log_link"))
summary(alpha_hat)

对于正态分布，以下计算确实给了我对 sd 的估计。

x <- rnorm(n = 10000, mean = 160, sd = 6)
theta_1 <- maxlogL(x = x, dist = 'dnorm', control = list(trace = 1),link = list(over = "sd", fun = "log_link"),
                 fixed = list(mean = 160))
summary(theta_1)

有人可以指出第一段代码中的错误吗？

4

1 回答 1

0

我不知道。我会变得懒惰，并以我更熟悉的不同方式来做：

library(bbmle)
m <- mle2(x~dbeta(shape1=exp(loga),shape2=1), 
        data=data.frame(x), start=list(loga=0))
## estimate (back-transformed)
exp(coef(m)) ## 0.6731152
## profile confidence interval (back-transformed)
exp(confint(m))
##     2.5 %    97.5 %
## 0.6322529 0.7157005

设置一个简单的引导函数...

bootfun <- function() {
    newx <- sample(x,size=length(x),replace=TRUE)
    newm <- update(m, data=data.frame(x=newx))
    return(coef(newm))
}
set.seed(101)
bootsamp <- replicate(500, bootfun())
exp(quantile(bootsamp, c(0.025, 0.975)))
##      2.5%     97.5% 
## 0.6533478 0.7300200

事实上，对于这种情况，（非常快的）Wald 置信区间可能很好......

exp(confint(m,method="quad"))
##     2.5 %    97.5 % 
## 0.6462591 0.7315456
于 2020-05-22T19:03:50.020 回答