如何优化包含 R 或 Python 中的积分的函数？

Question

第一次在这里提出问题，但绝对受益于浏览其他问题的所有惊人答案。

所以，问题是我正在尝试通过 R 中的 MLE 优化函数（或者 Python 也可以），并且函数的一部分是包含我正在估计的参数的积分。函数如下所示：

The log-likelihood function

其中，F(w) 是对数正态分布的 CDF，其均值和 sigma 未知，lambda 和 eta 也是如此。所有其他参数都是已知的，包括 w*。我尝试使用 R 中的 mle 命令估计此函数。但由于不可分割的部分，它无法通过。谁可以帮我这个事？我的代码如下。非常感谢！

Definitions of f(w) and F(w)

我的代码不起作用：

  func <- function(mu, sigma_w){
F_w = (-.5*(((log(w_star)-mu)/sigma_w))^2)

}

ll <- function(lambda,sigma_w,eta,mu) {
  R <- N*log(lambda) + N_u*log(1/sqrt(2*pi)*exp(integrate(func,lower = -Inf, upper = w_star, rel.tol=1e-5)$value)) - lambda*1/sqrt(2*pi)*exp(integrate(func,lower = -Inf, upper = w_star, rel.tol=1e-5)$value)*sum(T_ui) + N_u*log(eta) + sum(log(1/(sigma_w*w_i)*1/sqrt(2*pi)*exp(-.5*((log(w_i-mu)/sigma_w)^2)))) - N*log(eta+lambda*1/sqrt(2*pi)*exp(integrate(integrate(func,lower = -Inf, upper = w_star, rel.tol=1e-5)$value)))
  return(R)                           
}

fit <-mle(ll,start = list(lambda =0, sigma_w = 0, eta =0, mu = 0))

错误信息：

Error in f(x, ...) : argument "sigma_w" is missing, with no default 

Answer 1

这是在 R 中使用闭包（返回函数的函数）的一个很好的例子。它们正好为 MLE 服务，用你的值试试这个代码：

library(stats4)

ll <- function(lambda, eta, sigma_w, mu, w_start, w_i, N, N_u )
{
  # fix the others params
  w_start <- w_start
  mu <- mu
  T_ui <- T_ui #your sample
  w_i <- w_i #another sample
  sigma_w <- sigma_w
  N <- N
  N_u <- N_u 
  #return another funtion
  function(lambda, eta)
  {
    F_star <- log( plnorm( q = w_star,
                           meanlog = mu, sdlog = sigma_w) ) 
    ll.scalar <- N * lambda + N_u*log( F_star) - lambda * F_star * sum(T_ui) +
      N_u* log(eta) + sum( log( dlnorm(x = w_i))) - N * log(eta + F_star)

    return(-ll.scalar) # negative for use in stats4::mle
  }
}

# your values example
w_star <- 1
mu <- 1
T_ui <- c(1,2,3) #your sample
w_i <- c(1,2,3)/10 #another sample
sigma_w <- 1
N <- 1 
N_u <- 1 
# You build the function with the closure
myLL <- ll(lambda , eta, sigma_w , mu, w_star , w_i , N, N_u  )


# now optimitation 

result <- mle(myLL, start= list(lambda=0,  eta=0))
result

```