更新：最小限の例を下にスクロールします（同じエラーを再現しませんが、結果は異なります）
実際の質問の前に注意事項があります。私は R にまったく慣れていないので、理解できることを願っています。Google と R のデバッグ機能を使用して問題を突き止めようとしました (ただし、まだあまり慣れていません)。
私はとりわけ読んだことがあります：
サブセット化中に %in% と == をいつ使用するかについての質問
Hadley Wickham の非常に優れたブログにたどり着きました。ここでは、R のデバッグ機能について詳しく調べます。彼の記事: 非標準評価
私の問題
私は、とりわけ社債に関するデータを含む dataと呼ばれるデータフレームを持っています。
最適関数のラッパーである bbmle パッケージの mle2 関数を使用して、用語構造のモデルを推定したいと考えています。以下のコードを参照してください。正規分布の残差の仮定が正当化されるかどうかは疑問の余地がありますが、次のとおりです。
アプローチ 1: 「手動で」サブセット化しても問題なく動作します。
require("bbmle")
lifeBBB <- data$life[data$Rating == "BBB"]
yieldBBB <- data$yield[data$Rating == "BBB"]
LL <- function(b0, b1, b2, lambda, mu, sigma){
Score= yieldBBB - (b0+b1*((1-exp(-lambda*lifeBBB))/(lambda*lifeBBB))+ b2*((1-exp(-lambda*lifeBBB))/(lambda*lifeBBB)-exp(-lambda*lifeBBB)))
Score = suppressWarnings(dnorm(Score, mu, sigma, log=TRUE))
-sum(Score)
}
fit <- mle2(LL, start = list(b0 = 108, b1 = -85, b2=168, lambda= 0.12, mu = -10, sigma=70),control=list(maxit = 5000))'
Call:
mle2(minuslogl = LL, start = list(b0 = 108, b1 = -85, b2 = 168,
lambda = 0.12, mu = -10, sigma = 70), control = list(maxit = 5000))
Coefficients:
b0 b1 b2 lambda mu
110.2355968 -82.7010072 167.5960478 0.1410541 -7.7644032
sigma
69.4846302
Log-likelihood: -1018.8
アプローチ 2: mle2 関数内からサブセットを呼び出そうとすると、エラーが発生します。
LL2 <- function(b0, b1, b2, lmbd, mu, sigma){
Score= yield - (b0+b1*((1-exp(-lmbd*life))/(lmbd*life))+ b2*((1-exp(-lmbd*life))/(lmbd*life)-exp(-lmbd*life)))
Score = suppressWarnings(dnorm(Score, mu, sigma, log=TRUE))
-sum(Score)
}
fit1 <- mle2(minuslogl=LL2, start = list(b0 = 108, b1 = -85, b2=168, lmbd= 0.12, mu = -10, sigma=70),method="BFGS", data=data, subset= Rating=="BBB", control=list(maxit=5000))
Error in optim(par = c(108, -85, 168, 0.12, -10, 70), fn = function (p) :
initial value in 'vmmin' is not finite
したがって、アプローチ 1 は機能し、アプローチ 2 は、評価ごとに新しい変数を導入する必要がなく、まったく同じことを行うより便利な方法であるため、コードに何か問題があるに違いないと考えました。私が言ったように、私はデバッグ方法を使用しようとしましたが、R のマシンルームの奥深くで、まったく別のエラーが発生しました。
生成されたエラーについて、 R メーリング リストからこれを見つけました。ここでの問題は、dbinom 分布のサイズ パラメータとして整数ではなく分数が指定されていたことです。しかし、これが私のコードの問題であることがわかりません。
前もって感謝します
B.ローア
リクエストに応じて、模範的なデータを使用した最小限の例を考え出そうとしました。最小限の例ではエラーは発生しませんが、2 つの異なる推定値が生成されます。私の理解では、両方の関数が同じことをするはずなので、これは独特です
#### minmial example
require("bbmle")
## approach 1:
set.seed(23456)
Rating <- c(rep("A",38),rep("BBB",39) )
yield <- rnorm(77,0.5,15)
life <- runif(77,1,20)
exdata <- data.frame(Rating,yield,life)
lifeBBB <- exdata$life[exdata$Rating == "BBB"]
yieldBBB <- exdata$yield[exdata$Rating == "BBB"]
LL <- function(b0, b1, b2, lambda, mu, sigma){
Score= yieldBBB - (b0+b1*((1-exp(-lambda*lifeBBB))/(lambda*lifeBBB))+ b2*((1-exp(-lambda*lifeBBB))/(lambda*lifeBBB)-exp(-lambda*lifeBBB)))
Score = suppressWarnings(dnorm(Score, mu, sigma, log=TRUE))
-sum(Score)
}
fit <- mle2(LL, start = list(b0 = 100, b1 = -80, b2=160, lambda= 0.12, mu = 1, sigma=14),method ="BFGS",control=list(maxit = 5000))
### approach 2:
LL2 <- function(b0, b1, b2, lmbd, mu, sigma){
Score= exdata$yield - (b0+b1*((1-exp(-lmbd*exdata$life))/(lmbd*exdata$life))+ b2*((1-exp(-lmbd*exdata$life))/(lmbd*exdata$life)-
exp(-lmbd*exdata$life)))
Score = suppressWarnings(dnorm(Score, mu, sigma, log=TRUE)) ## assumption is that residuals are normally distributed
-sum(Score)
}
fit1 <- mle2(minuslogl=LL2, start = list(b0 = 100, b1 = -80, b2=160, lmbd= 0.12, mu = 1, sigma=14),method="BFGS", data=exdata, subset= Rating=="BBB", control=list(maxit=5000))
Call(fit):
Call:
mle2(minuslogl = LL, start = list(b0 = 100, b1 = -80, b2 = 160,
lambda = 0.12, mu = 1, sigma = 14), method = "BFGS", control = list(maxit = 5000))
Coefficients:
b0 b1 b2 lambda mu sigma
94.34166416 -85.35582080 159.76952349 -0.00283408 -4.65833584 12.94283927
Log-likelihood: -155.2
Call (fit1):
all:
mle2(minuslogl = LL2, start = list(b0 = 100, b1 = -80, b2 = 160,
lmbd = 0.12, mu = 1, sigma = 14), method = "BFGS", data = exdata,
subset = Rating == "BBB", control = list(maxit = 5000))
Coefficients:
b0 b1 b2 lmbd mu sigma
9.269496e+01 -8.600709e+01 1.586543e+02 -1.186371e-04 -6.305040e+00 1.458118e+01
Log-likelihood: -315.6