R優化多個參數

Question

我正在使用R optim（）函數來估計優化用戶定義函數的參數集，如下所示。但是的Optim（）出放是：

錯誤的Optim（PSTART，llAgedepfn，方法= 「L-BFGS-B」，上部=起來，低級= LO）： L-BFGS-B需要有限值的 'FN'

請幫助。完整的腳本如下所示：

dataM<-cbind(c(1.91,0.29,0.08,0.02,0.01,0.28,0.45,0.36,0.42,0.17,0.16,0.06,0.17,0.17,0.12), 
       c(0.27,4.54,0.59,0.05,0.04,0.13,0.48,0.68,0.66,0.18,0.11,0.06,0.08,0.08,0.08), 
       c(0.07,0.57,4.48,0.48,0.02,0.05,0.09,0.43,0.78,0.52,0.17,0.10,0.05,0.05,0.14), 
       c(0.02,0.04,0.44,4.34,0.36,0.09,0.07,0.11,0.41,0.77,0.43,0.10,0.03,0.04,0.14), 
       c(0.01,0.04,0.01,0.36,2.20,0.46,0.19,0.15,0.19,0.34,0.62,0.30,0.09,0.03,0.22), 
       c(0.22,0.11,0.05,0.09,0.45,0.91,0.61,0.43,0.37,0.26,0.41,0.63,0.29,0.16,0.15), 
       c(0.31,0.35,0.07,0.05,0.16,0.54,0.81,0.59,0.48,0.36,0.33,0.43,0.47,0.26,0.20), 
       c(0.22,0.45,0.29,0.08,0.11,0.34,0.53,0.85,0.71,0.39,0.27,0.26,0.26,0.28,0.38), 
       c(0.22,0.36,0.44,0.26,0.12,0.24,0.36,0.59,0.91,0.61,0.35,0.28,0.20,0.22,0.29), 
       c(0.09,0.10,0.30,0.49,0.22,0.17,0.28,0.33,0.62,0.80,0.52,0.29,0.20,0.11,0.46), 
       c(0.10,0.07,0.12,0.32,0.48,0.32,0.30,0.27,0.42,0.61,0.78,0.47,0.33,0.23,0.49), 
       c(0.04,0.04,0.06,0.08,0.24,0.53,0.41,0.28,0.36,0.36,0.50,0.67,0.51,0.19,0.47), 
       c(0.10,0.05,0.04,0.02,0.07,0.23,0.43,0.26,0.23,0.23,0.33,0.48,0.75,0.51,0.49), 
       c(0.05,0.04,0.03,0.05,0.02,0.10,0.19,0.22,0.21,0.10,0.18,0.14,0.40,0.79,0.82), 
       c(0.03,0.02,0.03,0.03,0.06,0.04,0.06,0.12,0.11,0.18,0.16,0.14,0.16,0.34,1.26) 
) 

NormCM <- dataM/eigen(CMWkday)$values[1] #Normalizing the contact mtrix - divide by the largest eigen value 

w <- c(495,528,548,603,617,634,720,801,957,937,798,755,795,1016,2469) 

g2 <- c(770,622,726,559,410,547,564,472,399,397,340,308,337,91,84) 

h2 <- c(269,426,556,430,271,284,303,207,194,181,126,106,74,24,23) 

z2 <- h2/g2 

g1 <- c(774,527,665,508,459,539,543,492,402,412,365,342,213,146,152) 

h1 <- c(56,31,84,173,103,85,123,70,71,80,55,25,18,12,26) 
z1 <- h1/g1 

#### Normal loglikelihood ######### 

llnormfn <- function(q) { 

    tol <- 1e-9 
    final.size.start <- 0.8 
    zeta <- rep(final.size.start, nrow(NormCM)) 
    last.zeta <- rep(0, nrow(NormCM)) 
    first.run <- T 
    current.diff <- tol+1 
    loglik <- 0 

    while (current.diff > tol) { 

    zeta <- 1-exp(-(q*(zeta%*%NormCM))) 
    current.diff <- sum(abs(last.zeta-zeta)) 
    last.zeta <-zeta 

    } 
    mu <- c(zeta) 

    zigma <- z1*(1-z1)/g1 + (z1+mu)*(1-(z1+mu))/g2 

    logliknorm <- -sum((((z2-z1)-mu)**2)/2*zigma + 0.5*log(2*pi*zigma)) 

    return(logliknorm) 

} 

pstart <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1) 
up <- c(5,5,5,5,5,5,5,5,5,5,5,5,5,5,5) 
lo <- c(0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1) 
estm <- optim(pstart, llnormfn, method = "L-BFGS-B", upper = up, lower = lo) 

Answer 1

您llnormfn不返回範圍內的所有參數的值的有限值。例如在上限：

> llnormfn(up) 
[1] NaN 
Warning message: 
In log(2 * pi * zigma) : NaNs produced 

因爲zigma在這裏必須小於零。

如果限制的範圍內位，你可以最終找到一個地方，它的工作...

> llnormfn(up-2) 
[1] NaN 
Warning message: 
In log(2 * pi * zigma) : NaNs produced 
> llnormfn(up-3) 
[1] 42.96818 

讓我們來看看它工作在較低的範圍內：

> llnormfn(lo) 
[1] 41.92578 

看起來罰款。因此，無論您是否將函數的上限設置在函數的計算有效範圍之外，或者您的函數或其他函數中都存在錯誤。

如果您運行與優化降低上限你得到收斂：

> estm <- optim(pstart, llnormfn, method = "L-BFGS-B", upper = up-3, lower = lo) 
> estm 
$par 
[1] 1.9042672 1.0891264 0.9916916 0.6208685 1.2413983 1.4822433 1.1243878 
[8] 1.5224263 1.3686933 1.4876350 1.6231518 2.0000000 2.0000000 2.0000000 
[15] 2.0000000 

$value 
[1] 38.32182 

$counts 
function gradient 
     23  23 

$convergence 
[1] 0 

$message 
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH" 

雖然你可能會注意到其中的一些參數是在的上限值（2.0），這是一個警鐘。

檢查功能很好地處理其輸入值 - 試着修改所有，但一和密謀如何llnormfn的行爲，同時改變之一。我只是看了一眼，功能看起來並不平坦，有很多不連續性，所以我認爲BFGS是一種優化的好方法。

例如改變在0.1和2之間的第五個參數：

> s = seq(0.1,2,len=300) 
> ss = sapply(1:length(s),function(i){ll=lo;ll[5]=s[i];llnormfn(ll)}) 
> plot(s,ss) 

給出：