R推定パラメータ

Question

私は統計パッケージから機能optimを使用していたが、あるR.

における最尤によって二項分布から推定パラメータNとPをしようとしていますエラー。

私のコードである：OPTIMで

xi = rbinom(100, 20, 0.5) # Sample 
n = length(xi) # Sample size 

# Log-Likelihood 
lnlike <- function(theta){ 
log(prod(choose(theta[1],xi))) + sum(xi*log(theta[2])) + 
(n*theta[1] - sum(xi))*log(1-theta[2]) 
} 

# Optimizing 
optim(theta <- c(10,.3), lnlike, hessian=TRUE) 

エラー（シータ< - C（10、0.3）、lnlike、ヘッセ= TRUE）： 関数は初期パラメータ

で評価することができません

誰でもこれを行いましたか？どの機能が使われましたか？

Answer 1

TL; DRあなたが応答変数は、理論上の最大値である二項N（より大きい場合はゼロの可能性（したがって、負の無限対数尤度）を取得するつもりですレスポンス）。実際の問題では、Nは既知のものとして扱われ、確率は推定されます。 Nを見積もりたいのであれば、（1）サンプルの最大値以上になるように制約する必要があります。 （2）離散的でなければならないパラメータに対して最適化するために特別なことをする（これは高度な/トリッキーな問題である）。

最初の部分は問題を特定するためのデバッグ戦略を示し、2番目はNとpを同時に最適化する戦略を示しています。

セットアップ：

set.seed(101) 
n <- 100 
xi <- rbinom(n, size=20, prob=0.5) # Sample 

対数尤度関数：

lnlike <- function(theta){ 
    log(prod(choose(theta[1],xi))) + sum(xi*log(theta[2])) + 
     (n*theta[1] - sum(xi))*log(1-theta[2]) 
} 

はこれを打破するのをしてみましょう。

theta <- c(10,0.3) ## starting values 
lnlike(c(10,0.3)) ## -Inf 

OK、対数尤度は、開始値で-Infです。驚くことではないがoptim()はそれで動作しません。

用語を解説しましょう。

log(prod(choose(theta[1],xi))) ## -Inf 

OK、私たちはすでに最初の言葉に困っています。

prod(choose(theta[1],xi)) ## 0 

製品はゼロです...なぜですか？

choose(theta[1],xi) 
## [1] 120 210 10 0 0 10 120 210 0 0 45 210 1 0 

ゼロがたくさんあります。どうして？ xiの値は何か問題がありますか？

## [1] 7 6 9 12 11 9 7 6 

アハ！私たちは、あなたが本当にこれをしたい場合は、あなたがnの妥当な値上ブルートフォース列挙によってそれを行うことができます12

badvals <- (choose(theta[1],xi)==0) 
all(badvals==(xi>10)) ## TRUE 

で7、6、9 ...しかし、トラブルにするためにOKです。.. 。

## likelihood function 
llik2 <- function(p,n) { 
    -sum(dbinom(xi,prob=p,size=n,log=TRUE)) 
} 
## possible N values (15 to 50) 
nvec <- max(xi):50 
Lvec <- numeric(length(nvec)) 
for (i in 1:length(nvec)) { 
    ## optim() wants method="Brent"/lower/upper for 1-D optimization 
    Lvec[i] <- optim(par=0.5,fn=llik2,n=nvec[i],method="Brent", 
        lower=0.001,upper=0.999)$val 
} 
nvec[which.min(Lvec)] ## 20 
par(las=1,bty="l") 
plot(nvec,Lvec,type="b") 

Answer 2

なぜトラブルになるのですか？

lnlike(c(10, 0.3))を入力すると、-Infとなります。あなたのエラーメッセージは、optimではなく、lnlikeです。

多くの場合、nが知られており、pしか推定する必要はありません。この状況では、モーメント推定器または最尤推定器のいずれかが閉鎖された形式であり、数値最適化は必要ない。ですから、nと推測するのは本当に変です。

推定する場合は、制約があることに注意する必要があります。

range(xi) ## 5 15 

あなたの観測範囲[5, 15]を持っているので、それはn >= 15ことが必要とされるチェック。どのようにして初期値10を渡すことができますか？ nの検索方向は大きい開始値から、max(xi)に達するまで徐々に下方に検索する必要があります。nの初期値として30を試してみてください。

さらに、現在の方法でlnlikeを定義する必要はありません。次の操作を行います（それが最大化を行うことができますが）

lnlike <- function(theta, x) -sum(dbinom(x, size = theta[1], prob = theta[2], log = TRUE)) 

• optimは、多くの場合、最小化のために使用されています。負の対数尤度を得るために関数にマイナス記号を入れました。このようにして、lnlike w.r.tを最小化しています。 theta。
• また、グローバル環境からの観測ではなく、の追加の引数として、観測値xiをlnlikeに渡す必要があります。

optimとナイーブ試し：

私のコメントでは、私はすでにoptimが連続変数のために使用されている間、私はnは整数でなければならないので、nが動作する推定するoptimを使用して信じていないことを言いました。これらのエラーと警告はあなたに納得させるでしょう。

optim(c(30,.3), fn = lnlike, x = xi, hessian = TRUE) 

Error in optim(c(30, 0.3), fn = lnlike, x = xi, hessian = TRUE) : 
    non-finite finite-difference value [1] 
In addition: There were 15 or more warnings (use warnings() to see the 
first 15 

> warnings() 
Warning messages: 
1: In dbinom(x, size = theta[1], prob = theta[2], log = TRUE) : NaNs produced 
2: In dbinom(x, size = theta[1], prob = theta[2], log = TRUE) : NaNs produced 
3: In dbinom(x, size = theta[1], prob = theta[2], log = TRUE) : NaNs produced 
4: In dbinom(x, size = theta[1], prob = theta[2], log = TRUE) : NaNs produced 
5: In dbinom(x, size = theta[1], prob = theta[2], log = TRUE) : NaNs produced 

解決策？

ベンはあなたに道を提供しました。 をnと見なす代わりに、手動でnのグリッド検索を行います。 nの各候補について、単変量最適化w.r.tを実行します。 p。 実際には数値最適化を行う必要はありません）このようにして、プロファイルの可能性はnになります。次に、グリッド上にnがあり、このプロファイルの可能性を最小限に抑えることができます。

ベンは完全な詳細を提供しており、私はそれを繰り返さないものとします。ニース（そして素早く）仕事、ベン！