有没有办法让已经定义的参数在 R 中使用 optim() 函数时显示为丢失?
Is there a way so that an already defined argument can appear as missing when using optim() function in R?
我正在尝试获取 Gumbel 分布的对数似然的最大似然估计量以进行生存分析(我这么说是为了让您不会对对数似然函数感到陌生,我认为它是正确的)。为了做到这一点,我必须通过使用 optim 函数来最大化负对数似然,我尝试这样做但控制台在 fn(par, ...) 中给我一个错误:参数 "b" 是缺失,没有默认值。
我也尝试过以与此 link 的答案类似的方式进行操作:Solve for maximum likelihood with two parameters under constraints,但控制台游戏如下:
optim(c(1, 1), function(x) log_lhood(x[1], x[2], d = lung$status, 错误:
objective optim 中的函数求值为长度 0 而不是 1。
log_lhood <- function(m,b,d,t){
sum<-0
for (i in 1:length(lung)){
if (d[i] == 1){
sum<- sum - log(1-exp(-exp(-(t[i]-m)/b)))
} else {
sum<- sum - log((1/b)*exp(-(t[i]-m/b+exp(-(t[i]-m/b)))))
}
}
}
#a,b parameter optimization
optim(c(0,0), fn = log_lhood, d = lung$status, t = KM_fit$time) #1st way
optim(c(1, 1), function(x) log_lhood(x[1], x[2],d=lung$status,t=KM_fit_test$time)) #2nd way as in the link
这里有一些问题。该函数的第一个参数应该是一个参数向量。另外你需要 nrow(lung)
而不是 length(lung)
,最好使用 length(d)
。另外你不应该在这里使用循环,它非常低效,使用 ifelse()
(在 R 中我们总是尝试向量化一切)。您还需要检查是否可以为所有参数值计算对数似然(例如 b=0)。你也忘了return(sum)
。另外 sum
是一个你不应该屏蔽的有用函数。
这会运行。
library(reprex)
lung <- data.frame(status=c(0,0,1,1))
KM_fit <- data.frame(time=c(0,1,2,3))
log_lhood <- function(x,d,t){
m <- x[1]
b <- x[2]
sum <-0
for (i in 1:nrow(lung)){
if (d[i] == 1){
sum <- sum - log(1-exp(-exp(-(t[i]-m)/b)))
} else {
sum <- sum - log((1/b)*exp(-(t[i]-m/b+exp(-(t[i]-m/b)))))
}
}
return(sum)
}
#a,b parameter optimization
optim(par=c(0,1), fn = log_lhood, d = lung$status, t = KM_fit$time)
$par
[1] 1.661373 1.811780
$value
[1] 5.318068
$counts
function gradient
63 NA
$convergence
[1] 0
$message
NULL
您可以像这样重写您的函数。
log_lhood <- function(x,d,t){
m <- x[1]
b <- x[2]
s <- ifelse(d==1,
-log(1-exp(-exp(-(t-m)/b))),
-log((1/b)*exp(-(t-m/b+exp(-(t-m/b)))))
)
return(sum(s, na.rm=TRUE))
}
我正在尝试获取 Gumbel 分布的对数似然的最大似然估计量以进行生存分析(我这么说是为了让您不会对对数似然函数感到陌生,我认为它是正确的)。为了做到这一点,我必须通过使用 optim 函数来最大化负对数似然,我尝试这样做但控制台在 fn(par, ...) 中给我一个错误:参数 "b" 是缺失,没有默认值。
我也尝试过以与此 link 的答案类似的方式进行操作:Solve for maximum likelihood with two parameters under constraints,但控制台游戏如下: optim(c(1, 1), function(x) log_lhood(x[1], x[2], d = lung$status, 错误: objective optim 中的函数求值为长度 0 而不是 1。
log_lhood <- function(m,b,d,t){
sum<-0
for (i in 1:length(lung)){
if (d[i] == 1){
sum<- sum - log(1-exp(-exp(-(t[i]-m)/b)))
} else {
sum<- sum - log((1/b)*exp(-(t[i]-m/b+exp(-(t[i]-m/b)))))
}
}
}
#a,b parameter optimization
optim(c(0,0), fn = log_lhood, d = lung$status, t = KM_fit$time) #1st way
optim(c(1, 1), function(x) log_lhood(x[1], x[2],d=lung$status,t=KM_fit_test$time)) #2nd way as in the link
这里有一些问题。该函数的第一个参数应该是一个参数向量。另外你需要 nrow(lung)
而不是 length(lung)
,最好使用 length(d)
。另外你不应该在这里使用循环,它非常低效,使用 ifelse()
(在 R 中我们总是尝试向量化一切)。您还需要检查是否可以为所有参数值计算对数似然(例如 b=0)。你也忘了return(sum)
。另外 sum
是一个你不应该屏蔽的有用函数。
这会运行。
library(reprex)
lung <- data.frame(status=c(0,0,1,1))
KM_fit <- data.frame(time=c(0,1,2,3))
log_lhood <- function(x,d,t){
m <- x[1]
b <- x[2]
sum <-0
for (i in 1:nrow(lung)){
if (d[i] == 1){
sum <- sum - log(1-exp(-exp(-(t[i]-m)/b)))
} else {
sum <- sum - log((1/b)*exp(-(t[i]-m/b+exp(-(t[i]-m/b)))))
}
}
return(sum)
}
#a,b parameter optimization
optim(par=c(0,1), fn = log_lhood, d = lung$status, t = KM_fit$time)
$par
[1] 1.661373 1.811780
$value
[1] 5.318068
$counts
function gradient
63 NA
$convergence
[1] 0
$message
NULL
您可以像这样重写您的函数。
log_lhood <- function(x,d,t){
m <- x[1]
b <- x[2]
s <- ifelse(d==1,
-log(1-exp(-exp(-(t-m)/b))),
-log((1/b)*exp(-(t-m/b+exp(-(t-m/b)))))
)
return(sum(s, na.rm=TRUE))
}