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 

正常に見えます。だからあなたの関数の計算上有効な範囲外にその上限を設定したか、またはllnormfn関数、あるいはその両方にバグがあります。

あなたは減少し、上側あなたが収束得るかバインドして最適化を実行しない場合：あなたはそれらのパラメータの一部は警鐘であるの上限値（2.0）です気づくかもしれませんが

> 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" 

を。

関数が入力値に対して賢明に動作するかどうかを確認してください。オールと1を固定して、どのように変化させるかをプロットしてください。llnormfn私はちょうど見た目が素早く、関数は滑らかに見えず、多くの不連続点があるので、BFGSは最適化のための良い方法であるとは思わない。

例えば、0.1と2の間の5番目のパラメータを変化させる：

> 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) 

を与える：