你能帮我解决R中optim函数的错误吗?

Question

我是r的新手，我想用最大似然估计做一些参数估计。

这是我的尝试:

数据

my_data = c(0.1,0.2,1,1,1,1,1,2,3,6,7,11,12,18,18,18,18,18,21,32,36,40,
45,45,47,50,55,60,63,63,67,67,67,67,72,75,79,82,82,83,
84,84,84,85,85,85,85,85,86,86)

和

lx <- function(p,x){
l <- p[1]
b <- p[2]
a <- p[3]
n <- length(x)
lnL <- n*log(l)+n*log(b)+n*log(a)+(b-1)*sum(log(x))+(a-1)*sum(log(1+l*x^b))+n-sum(1+l*x^b)
return(-lnL)
}

注:l为λ， b为β， a为α。

这是optim函数

optim(p=c(1,1,1),fn = lx, method = "L-BFGS-B",
lower = c(0.0001, 0.0001, 0.0001),
control = list(), hessian = FALSE, x = my_data)

在我运行这段代码之后，我得到一个错误消息:

Error in optim(p = c(1, 1, 1), fn = lx, method = "L-BFGS-B", lower = c(1e-04,  : 
objective function in optim evaluates to length 50 not 1

我的代码有什么问题?你能帮我修一下吗?提前感谢!

Answer 1

用MASS::fitdistr代替对数似然。

#
# Power Generalized Weibull distribution
#
# x > 0, alpha, beta, lambda > 0
#
dpowergweibull <- function(x, alpha, beta, lambda){
f1 <- lambda * beta * alpha
f2 <- x^(beta - 1)
f3 <- (1 + lambda * x^beta)^(alpha - 1)
f4 <- exp(1 - (1 + lambda * x^beta)^alpha)
f1 * f2 * f3 * f4
}
ppowergweibull <- function(q, alpha, beta, lambda){
1 - exp(1 - (1 + lambda * q^beta)^alpha)
}
my_data <- c(0.1,0.2,1,1,1,1,1,2,3,6,7,11,12,18,18,18,18,18,21,32,36,40,
45,45,47,50,55,60,63,63,67,67,67,67,72,75,79,82,82,83,
84,84,84,85,85,85,85,85,86,86)
start_par <- list(alpha = 0.1, beta = 0.1, lambda = 0.1)
y1 <- MASS::fitdistr(my_data, dpowergweibull, start = start_par),
start_par2 <- list(shape = 1, rate = 1)
y2 <- MASS::fitdistr(my_data, "gamma", start = start_par2)
hist(my_data, freq = FALSE)
curve(dpowergweibull(x, y1$estimate[1], y1$estimate[2], y1$estimate[3]),
from = 0.1, to = 90, col = "red", add = TRUE)
curve(dgamma(x, y2$estimate[1], y2$estimate[2]),
from = 0.1, to = 90, col = "blue", add = TRUE)