如何找到 nlimb 优化的最佳起始值？答案

【问题标题】：How do I find best starting values for nlimb optimization?如何找到 nlimb 优化的最佳起始值？
【发布时间】：2021-11-09 21:42:57
【问题描述】：

我正在尝试使用附加的数据集运行以下代码。如何解决hessian矩阵求逆的错误？

library(stats4)
library(bbmle)
library(stats)
library(numDeriv)
library('bbmle')
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)
hist(x)
fEHLKUMW<-function(a,b,alpha,vartheta)
{
  -sum(log(  (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)))
))
}
EHLKUMW.result<-mle2(fEHLKUMW,hessian = NULL,start=list(a=0.01,b=0.01,alpha=.3,vartheta=0.01),optimizer="nlminb",lower=0)
summary(EHLKUMW.result)

我得到的错误是；

**
Warning messages:
1: In nlminb(start = start, objective = objectivefunction, hessian = NULL,  :
  NA/NaN function evaluation
2: In mle2(fEHLKUML, hessian = NULL, start = list(a = 1, b = 0.4, c = 0.5,  :
  couldn't invert Hessian
**

【问题讨论】：

标签： r optimization parameters

【解决方案1】：

我认为这是一个非常开放的问题，因此我将介绍一些工具和方法，也许其他人可以发表评论等。

一、主要部分：

library(bbmle)
library(stats)
library(numDeriv)
library(bbmle)
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)
hist(x)
fEHLKUMW <- function(a,b,alpha,vartheta) {
  -sum(log(  (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)))
  ))
}

现在，我们当然可以像您一样运行它：

EHLKUMW.result <- mle2(
  fEHLKUMW,
  hessian = NULL,
  start = list(
    a = 0.01,
    b = 0.01,
    alpha = .3,
    vartheta = 0.01
  ),
  optimizer = "nlminb",
  lower = 0
)

但我们也可以在每个参数上使用分布来运行它，以始终获得新的输入：


EHLKUMW.result <- mle2(
  fEHLKUMW,
  hessian = NULL,
  start = 
    list(
      # a = 0.01,
      # a = rt(1, 10, ncp = 0.01),
      a = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10),
      # b = 0.01,
      b = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10),
      # alpha = .3,
      alpha = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.3, df = 10),
      # vartheta = 0.01
      vartheta = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10)
    ),    
  optimizer = "nlminb",
  lower = 0
)

我选择使用trundist 来获得t-分布，集中在您提供的值，较低的是 0 通过 a = -argument。如果您知道这些参数的上限是多少，可以使用b = -argument 来完成。

我认为最相关的输出是 logLik 和 coef。


library(tidyverse)

exec(fEHLKUMW, !!!list(
  a = 0.01,
  b = 0.01,
  alpha = .3,
  vartheta = 0.01
))
replicate(
  250,
  exec(fEHLKUMW, !!!list(
    # a = 0.01,
    # a = rt(1, 10, ncp = 0.01),
    a = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10),
    # b = 0.01,
    b = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10),
    # alpha = .3,
    alpha = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.3, df = 10),
    # vartheta = 0.01
    vartheta = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10)
  )))

我用它来微调参数的分布

以下是多次运行和它们的输出对比 logLik.


tibble(n = seq_len(100),
       output = map(n, ~mle2(
         fEHLKUMW,
         hessian = NULL,
         start = 
           list(
             # a = 0.01,
             # a = rt(1, 10, ncp = 0.01),
             a = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10),
             # b = 0.01,
             b = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10),
             # alpha = .3,
             alpha = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.3, df = 10),
             # vartheta = 0.01
             vartheta = truncdist::rtrunc(1, spec = "t", a = 0, ncp = 0.01, df = 10)
           ),    
         optimizer = "nlminb",
         lower = 0
       ))) ->
  outputs_df

这段代码打印效果很好

outputs_df %>% 
  mutate(coef = output %>% map(coef),
         logLik = output %>% map_dbl(logLik)) %>% 
  unnest_wider(coef) %>% 
  arrange(logLik) %>% 
  print(n=Inf)

【讨论】：