0

与 R 中的 optim 函数有关的问题

到目前为止,我有以下代码,需要知道输入我的 X 和 T 的值。X 是 10 个值的向量,T 是与均值和方差相关的 10*2 值的向量。我希望输出采用 alpha、mean1、mean2、var1 和 var2 的一个新值的格式。不确定如何正确获取输入数据。

我想在这个函数中运行 X 的所有值,但只运行 T 的第一行(10 个值),我不确定如何为 T 执行此操作。我对第二行有不同的函数。

R <-runif(10, 0, 1)
S <-1-R
T <-t(cbind(R,S))


X <- runif(10, 25, 35)


Data1 <- function(xy) { 
  alpha <- xy[1]
  mean1 <- xy[2]
  mean2 <- xy[3]
  var1 <- xy[4]
  var2 <- xy[5] 

  -sum(0.5*(((X)-mean1)/var1)^2+alpha*mean1+log(2.5*var1)+log(exp(-alpha*mean1)+exp(-alpha*mean2))*(T))
}
starting_values <- c(0.3, 28, 38, 4, 3)
optim(starting_values, Data1, lower=c(0, 0, 0, 0, 0), method='L-BFGS-B')

还收到以下错误代码:

Error in optim(starting_values, Data1, lower = c(0, 0, 0, 0, 0), method = "L-BFGS-B") : 
  L-BFGS-B needs finite values of 'fn'

为任何帮助而欢呼。

编辑

包含的第二个功能

0.5*((y1-mean2)/var2)^2+alpha*mean2+log(2.5*var2)+ log(exp(-alpha*mean1)+exp(-alpha*mean2)))*T

好的,尽可能清楚地解释我想要做什么。上面原始帖子中的第一个函数一次获取 X 的所有 10 个值,并且应该获取第一行 T 数据(此处标记为 R),我不确定如何执行此操作。

上面详述的第二个函数应该再次连续获取所有 10 个 X 值,然后是来自 T 的第二行数据(下面标记为 S)

然后将所有这些加在一起。因此估计了五个未知参数。

       [,1]      [,2]      [,3]       [,4]      [,5]        [,6]      [,7]      [,8]      [,9]     [,10]
R 0.1477715 0.3055021 0.2963543 0.04149945 0.8342484 0.996865333 0.1592568 0.4623762 0.8000778 0.6979342
S 0.8522285 0.6944979 0.7036457 0.95850055 0.1657516 0.003134667 0.8407432 0.5376238 0.1999222 0.3020658

编辑2

即使运行相同的种子,我也没有得到与本相同的值。我已经检查过我是否安装了所有软件包,并且看起来我已经安装了。我没有得到相同的最终答案,我也无法调用 opt2$par 的单个项目。我将提供前几行和最后几行,而不是提供大量输出。

0.3 28 38 4 3 -74.97014 -120.7212 
Loading required package: BB
Loading required package: quadprog
Loading required package: ucminf
Loading required package: Rcgmin
Loading required package: Rvmmin

Attaching package: ‘Rvmmin’

The following object(s) are masked from ‘package:optimx’:

    optansout

Loading required package: minqa
Loading required package: Rcpp
0.3 28 38 4 3 -74.97014 -120.7212 
0.9501 28 38 4 3 -176.3368 -265.9074 
1.9001 28 38 4 3 -324.7782 -478.4652 
0.9501 28.95 38 4 3 -179.9994 -260.8711 
0.9501 28 38.95 4 3 -176.3366 -283.0445 
0.9501 28 38 4.95 3 -176.7836 -265.9074 
0.9501 28 38 4 3.95 -176.3368 -254.6188 

.....................

16.32409 27.86113 38.54337 3.940143 2.504167 -2566.194 -3826.233 
16.32409 27.86113 38.54337 3.940044 2.504167 -2566.194 -3826.233 
16.32409 27.86113 38.54337 3.940093 2.504199 -2566.194 -3826.232 
16.32409 27.86113 38.54337 3.940093 2.504136 -2566.194 -3826.234 
> opt2$par
$par
[1] 16.324085 27.861134 38.543373  3.940093  2.504167

> opt2$par["mean1"]
$<NA>
NULL
4

1 回答 1

7

第一次破解:我在上面使用了你的代码。我set.seed(101)在开头添加了可重复性。

为了清楚起见,稍微重新格式化了函数,但没有改变任何重要的内容,并添加了一条cat()用于调试目的的语句:

Data1 <- function(xy) {
    alpha <- xy[1]; mean1 <- xy[2]; mean2 <- xy[3]
    var1 <- xy[4]; var2 <- xy[5]

    r1 <- -sum(0.5*((X-mean1)/var1)^2+
           alpha*mean1+
           log(2.5*var1)+
           log(exp(-alpha*mean1)+
               exp(-alpha*mean2))*T[1,])
    r2 <- -sum(0.5*((X-mean2)/var2)^2+
           alpha*mean2+
           log(2.5*var2)+
           log(exp(-alpha*mean1)+exp(-alpha*mean2))*T[2,])

    cat(xy,r1,r2,"\n")
   r1+r2
}

一个略微压缩的版本,(1)利用with使功能更清洁;(2) 使用 R 的复制和向量回收能力

Data2 <- function(xy) {
    with(as.list(xy),
    {
       mmat <- rep(c(mean1,mean2),each=ncol(T))
       vmat <- rep(c(var1,var2),each=ncol(T))
       -sum(0.5*((X-mmat)/vmat)^2+
          alpha*mmat+
          log(2.5*vmat)+
          log(exp(-alpha*mean1)+exp(-alpha*mean2))*T)
    })
}

现在我们需要一个命名的起始值向量:

 starting_values <- c(alpha=0.3, mean1=28, mean2=38, var1=4, var2=3)

检查结果是否匹配:

 Data1(starting_values) ##  [1] -195.6913
 Data2(starting_values) ##  [1] -195.6913

这失败了(但向我们提供了有关它如何失败的信息):

 optim(par=starting_values, Data1, lower=rep(1e-4,5), method='L-BFGS-B',
     control=list(trace=6))

它产生大量输出,以:

##  21.29998 27.97361 37.98915 4.011199 3.001 -6014.225 
## 21.29998 27.97361 37.98915 4.011199 2.999 -6014.225 
## 85.29991 27.89318 37.95606 4.04533 3 Inf 
## Error in optim(par = starting_values, Data1, lower = rep(1e-04, 5), 
##    method = "L-BFGS-B",  : 
##     L-BFGS-B needs finite values of 'fn'

这至少告诉你哪里出了问题。我现在将尝试逐个评估您的表达式,以查看哪个位溢出。

正如聊天室的评论者(贾斯汀)所说,

您的第三个术语 log(exp(...) + exp(...)) 非常快地进入 -Inf,因为 alpha、mean1 和 mean2 是无界的。exp(-大数 * 大数) ~ 0

为了进一步调试,您可以:

  • 尝试重新安排对函数的评估以避免溢出
  • 设置一些参数的上限以避免溢出
  • 具有函数测试并在适当的情况下返回非常大的值而不是 Inf

不幸的是,L-BFGS-B它比其他一些优化器更脆弱,并且不允许非有限值。

接下来我尝试了包中的bobyqa优化器optimx,它允许边界处理非有限值(并且是一种无导数的方法,通常比基于导数的方法稍慢但更健壮):它似乎有效好的,虽然我不知道答案是否合理。

library(optimx)
opt2 <- optimx(par=starting_values, 
      Data1, lower=rep(1e-4,5), method='bobyqa')
opt3 <- optimx(par=starting_values, 
      Data2, lower=rep(1e-4,5), method='bobyqa')

看起来不错(只要这是一个明智的答案,我不知道)。

> opt2$par
$par
    alpha     mean1     mean2      var1      var2 
16.330752 27.815324 38.497483  3.894179  2.447219 

> opt3$par
$par
    alpha     mean1     mean2      var1      var2 
16.330900 27.820813 38.491290  3.887975  2.456052 

请注意,答案略有不同(在 var2 的情况下大约为 0.5%),这表明拟合可能稍微不稳定/表面可能非常平坦。(Data1和 Data2 应该给出相同的答案,并为起始值这样做,但我猜操作的顺序使它们对某些输入给出非常不同的答案 - 或者我在某个地方搞砸了......)

要从此拟合中提取单个分量,例如mean1,使用向量索引:

opt3$par["mean1"]  ## 27.820813
于 2012-08-23T15:17:08.153 回答