カスタムメイドの負の指数確率密度関数を使用して、いくつかのフィールド データを適合させたいと考えています。(動機 - 最終的には、 Bullock Shea and Skarpaas 2006の表 3 の多くの分布にフィールド データを当てはめたいと考えています)。
最初にdnegexp
、この投稿に従って関数を定義しました:
mle2 数式呼び出しのカスタム密度関数定義でエラーが発生しました
dnegexp <- function(x, mya, myb, log=FALSE){
logne <- log(mya)-myb*x
if(log) return(logne) else return(exp(logne))
}
次にrnegexp
、任意のパラメーターを使用して、この分布からデータセットを生成する関数を作成しました。(この機能の軽薄さをお許しください。この記事のトピックではありませんが、より良いものにするためのアイデアは大歓迎です...)
rnegexp <- function(n, thisa, thisb){
rnums <- array(data=NA,dim=n)
counter=0
while(counter<=n){
candidate <- runif(1,0,100)
if(runif(1,0,1)<dnegexp(candidate,thisa, thisb)) {
rnums[counter] <- candidate
counter <- counter+1
}
}
return(rnums)
}
これら 2 つの関数を使用して、既知のパラメーターを使用して負の指数分布に従うデータセットを作成して、ワークフローをテストします。
set.seed(501)
mynegexpvals <- rnegexp(100,0.08, 0.09)
hist(mynegexpvals, freq=FALSE)
mydists <- seq(0,100, by=1)
lines(mydists,dnegexp(mydists, 0.08, 0.09), col="blue", lwd=2)
パッケージから使用mle2
しbbmle
てパラメーターを見つけようとすると、開始値として正確な生成パラメーターを指定したにもかかわらず、無意味な値が返されます。
> library(bbmle)
> mle2(mynegexpvals ~ dnegexp(a,b), start = list(a=0.08, b=0.09), data=data.frame(mynegexpvals))
Call:
mle2(minuslogl = mynegexpvals ~ dnegexp(a, b), start = list(a = 0.08,
b = 0.09), data = data.frame(mynegexpvals))
Coefficients:
a b
2.421577e+12 -6.849330e+12
Log-likelihood: 6.148807e+15
検索空間を制限しようとすると、境界でパラメーター推定値が得られます。
> mle2(mynegexpvals ~ dnegexp(a,b), start = list(a=0.08, b=0.09), data=data.frame(mynegexpvals), method="L-BFGS-B", lower=c(a=0.04, b=0.0001), upper=c(a=1000, b=1000))
Call:
mle2(minuslogl = mynegexpvals ~ dnegexp(a, b), start = list(a = 0.08,
b = 0.09), method = "L-BFGS-B", data = data.frame(mynegexpvals),
lower = c(a = 0.04, b = 1e-04), upper = c(a = 1000, b = 1000))
Coefficients:
a b
1e+03 1e-04
Log-likelihood: 690.69
Warning message:
In mle2(mynegexpvals ~ dnegexp(a, b), start = list(a = 0.08, b = 0.09), :
some parameters are on the boundary: variance-covariance calculations based on Hessian may be unreliable
ここで大きな何かが欠けていますか?dnegexp 関数を間違って定義しましたか?