R中的优化问题

Question

与R中的优化函数有关的问题

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

我想在这个函数中运行X的所有值,但只有T的第一行(10个值),我不知道如何为T执行此操作.我对第2行有不同的功能.

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)并且我不知道如何执行此操作.

上面详述的第二个函数应该再次取X的所有10个值,然后从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

EDIT2

我没有获得与ben相同的值,即使运行相同的种子.我已经检查过我安装了所有软件包,看起来我会这么做.我没有得到相同的最终答案,我也无法调用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

Answer 1

第一个破解:我使用了上面的代码.我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(...))非常快,因为alpha,mean1和mean2是无界的.exp(-large number*large number)~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