如何找到下面函数中参数的 MLE 估计值?
How can I find MLE estimates for the parameters in the function below?
我有一个新函数,我想估计它的参数; a,b,alpha,vartheta 使用 MLE。我如何进行估算?
EMHL<-function(a,b,alpha,vartheta) {
(2*a*b*alpha*vartheta*
(x^(vartheta-1))* exp(-x^vartheta) *
((1-exp(-x^vartheta))^(a-1)) *
((1 - (1 -( 1 - exp(-x^vartheta))^a)^b)^(alpha-1)) *
(1- (1 - exp(-x^vartheta))^a )^(-b*(alpha+1))
}
对于给定的数据集
x<- c(1.1, 1.4, 1.3, 1.7,1.9, 1.8, 1.6, 2.2,
1.7, 2.7, 4.1, 1.8, 1.5, 1.2, 1.4, 3,
1.7, 2.3, 1.6, 2.0)
问题中显示的函数有语法错误,所以我们使用了最后注释中的那个。
如果该函数是一个密度函数并且您想最小化相应的负对数似然,那么尝试几个不同的起始值似乎会导致收敛。
(对于问题下方评论中给出的 x,list(a = 1, b = 1, alpha = 575, vartheta = 0.01)
似乎可以作为起始值。)
NLL <- function(par) -sum(log(EMHL(par[1], par[2], par[3], par[4])))
st <- list(a = 1, b = 1, alpha = 525, vartheta = .2)
res <- optim(st, NLL); res
给予:
$par
a b alpha vartheta
8.845296e-01 1.211526e+00 5.315759e+02 1.326975e-03
$value
[1] -10904.36
$counts
function gradient
327 NA
$convergence
[1] 0
$message
NULL
L-BFGS-B
将L-BFGS-B与上述函数一起使用会出现问题,但如果我们限制参数就可以获得答案。例如,这会在所有约束都处于活动状态时收敛,即解在可行域的边界上。请注意,从边界估计得出的检验和 p 值可能无效。
optim(c(1, 1, 1, 1), NLL, method = "L-BFGS-B", lower = 0.01, upper = 2)
其他密度
另一种可能是使用不同的分布。 descdist 制作的 Cullen & Frey 图表明伽马分布可能接近
library(fitdistrplus)
descdist(x)
(图后续)
或者我们可以试试 flexsurv 包中的广义伽玛 (dgengamma)。
library(bbmle)
library(flexsurv)
NLLgeng <- function(mu, sigma, Q) -sum(dgengamma(x, mu, sigma, Q, log = TRUE))
m <- mle2(NLLgeng, list(mu = 1, sigma = 1, Q = 1), optimizer = "nlminb")
summary(m)
library(fitdistrplus)
fit <- fitdist(x, "gengamma", start = list(mu = 1, sigma = 1, Q = 1))
summary(fit)
plot(fit)
备注
EMHL<-function(a,b,alpha,vartheta) {
2 * a * b * alpha * vartheta *
(x^(vartheta-1))* exp(-x^vartheta) *((1-exp(-x^vartheta))^(a-1)) *
((1 - (1 -( 1 - exp(-x^vartheta))^a)^b)^(alpha-1)) *
(1- (1 - exp(-x^vartheta))^a )^(-b*(alpha+1))
}
x<- c(1.1, 1.4, 1.3, 1.7,1.9, 1.8, 1.6, 2.2, 1.7, 2.7, 4.1, 1.8,
1.5, 1.2, 1.4, 3, 1.7, 2.3, 1.6, 2.0)
我有一个新函数,我想估计它的参数; a,b,alpha,vartheta 使用 MLE。我如何进行估算?
EMHL<-function(a,b,alpha,vartheta) {
(2*a*b*alpha*vartheta*
(x^(vartheta-1))* exp(-x^vartheta) *
((1-exp(-x^vartheta))^(a-1)) *
((1 - (1 -( 1 - exp(-x^vartheta))^a)^b)^(alpha-1)) *
(1- (1 - exp(-x^vartheta))^a )^(-b*(alpha+1))
}
对于给定的数据集
x<- c(1.1, 1.4, 1.3, 1.7,1.9, 1.8, 1.6, 2.2,
1.7, 2.7, 4.1, 1.8, 1.5, 1.2, 1.4, 3,
1.7, 2.3, 1.6, 2.0)
问题中显示的函数有语法错误,所以我们使用了最后注释中的那个。
如果该函数是一个密度函数并且您想最小化相应的负对数似然,那么尝试几个不同的起始值似乎会导致收敛。
(对于问题下方评论中给出的 x,list(a = 1, b = 1, alpha = 575, vartheta = 0.01)
似乎可以作为起始值。)
NLL <- function(par) -sum(log(EMHL(par[1], par[2], par[3], par[4])))
st <- list(a = 1, b = 1, alpha = 525, vartheta = .2)
res <- optim(st, NLL); res
给予:
$par
a b alpha vartheta
8.845296e-01 1.211526e+00 5.315759e+02 1.326975e-03
$value
[1] -10904.36
$counts
function gradient
327 NA
$convergence
[1] 0
$message
NULL
L-BFGS-B
将L-BFGS-B与上述函数一起使用会出现问题,但如果我们限制参数就可以获得答案。例如,这会在所有约束都处于活动状态时收敛,即解在可行域的边界上。请注意,从边界估计得出的检验和 p 值可能无效。
optim(c(1, 1, 1, 1), NLL, method = "L-BFGS-B", lower = 0.01, upper = 2)
其他密度
另一种可能是使用不同的分布。 descdist 制作的 Cullen & Frey 图表明伽马分布可能接近
library(fitdistrplus)
descdist(x)
(图后续)
或者我们可以试试 flexsurv 包中的广义伽玛 (dgengamma)。
library(bbmle)
library(flexsurv)
NLLgeng <- function(mu, sigma, Q) -sum(dgengamma(x, mu, sigma, Q, log = TRUE))
m <- mle2(NLLgeng, list(mu = 1, sigma = 1, Q = 1), optimizer = "nlminb")
summary(m)
library(fitdistrplus)
fit <- fitdist(x, "gengamma", start = list(mu = 1, sigma = 1, Q = 1))
summary(fit)
plot(fit)
备注
EMHL<-function(a,b,alpha,vartheta) {
2 * a * b * alpha * vartheta *
(x^(vartheta-1))* exp(-x^vartheta) *((1-exp(-x^vartheta))^(a-1)) *
((1 - (1 -( 1 - exp(-x^vartheta))^a)^b)^(alpha-1)) *
(1- (1 - exp(-x^vartheta))^a )^(-b*(alpha+1))
}
x<- c(1.1, 1.4, 1.3, 1.7,1.9, 1.8, 1.6, 2.2, 1.7, 2.7, 4.1, 1.8,
1.5, 1.2, 1.4, 3, 1.7, 2.3, 1.6, 2.0)