は、私がここにRで作業しています

Question

メトロポリス - ヘイスティングスアルゴリズム（MCMC）を使用して、プロビットモデルをシミュレートすることは私のコードです：私は問題があると思い

n <- 1000 
trueB0 <- 0 
trueB1 <- 0.1 
prior.mean0 <- 0 
prior.sd0 <- 10 
prior.mean1 <- 0 
prior.sd1 <- 10 

x <- rnorm(n, 0, 1) 
eta <- trueB0 + trueB1 * x 
P <- pnorm(eta) 
Y <- rbinom(n, 1, P) 

prior <- function(param){ 
    b0=param[1] 
    b1=param[2] 
    B0prior <- dnorm(b0, mean = prior.mean0, sd = prior.sd0, log = T) 
    B1prior <- dnorm(b1, mean = prior.mean1, sd = prior.sd1, log = T) 
    return(B0prior + B1prior) 
} 

likelihood <- function(param){ 
    b0=param[1] 
    b1=param[2] 
    eta.l <- b0 + b1 * x 
    P.l <- pnorm(eta.l) 
    singlelikelihoods = dbinom(Y, 1, P.l, log = T) 
    sumall=sum(singlelikelihoods) 
    return(sumall) 
} 

posterior <- function(param){ 
    return(prior(param)+likelihood(param)) 
} 

rho <- 0 
sd0 <- 1 
sd1 <- 1 
sigma <- matrix(c(sd0^2, sd0*sd1*rho, sd0*sd1*rho, sd1^2), ncol = 2) 

proposalfunction <- function(param){ 
    return(rmvnorm(1, mean = param, sigma)) 
} 

q <- function(m1, m2){ 
    return(dmvnorm(m1, mean = m2, sigma, log = T)) 
} 

run_metropolis_MCMC <- function(initialvalue, iterations){ 
    chain = array(dim = c(iterations+1,2)) 
    chain[1,] = initialvalue 
    for (i in 1:iterations){ 
    proposal = proposalfunction(chain[i, ]) 
    #bug here in the q evaluatoin part 
    alpha = exp(posterior(proposal) - posterior(chain[i, ]) + q(chain[i, ], proposal) - q(proposal, chain[i, ])) 
    if (runif(1) < min(1,alpha)) 
    { 
     chain[i+1,] = proposal 
    } 
    else 
    { 
     chain[i+1,] = chain[i,] 
    } 
    } 
    return(chain) 
} 

initialvalue <- c(0 , 0.1) 
iterations <- 10000 
MH = run_metropolis_MCMC(initialvalue, iterations) 

MH部のアルファを計算するときに、私がすべきことβ（β|β*）/ q（β* |β）（qは提案関数、βは現在の連鎖値、β*はβを用いてqによって生成された値）を得る。 しかし、私は実際にどのようにRコードでそれを得ることができるか分かりません。 （私はqが二変量正規分布密度でなければならないことを知っていますが、ベータとベータ*を関数に実装してqの評価を得る方法はわかりません）

Answer 1

私はあなたのコードを実行します。 q評価のバグはありません。あなたのq関数が正確にそれを与えるので、私は実際にあなたの質問が混乱しています：二変量法線（それはログ密度を与える）に対してq（鎖[i、]、proposal）= q（ベータ|ベータ*

一般的なカップルのあなたのコードのコメント：

• runif（1）<分（1、アルファ）はrunifと同じになります（1）<アルファ

• あなたの場合、それは良いことですコードは受け入れられた提案の数を数えるので、その割合を知ることができます。そのパーセンテージが小さい場合、提案は大きくなり、その逆もあります。あなたのforループの前にn.accept < -0を追加し、最初のif文でn.accept < - n.accept + 1を追加します。そして、リターン

  印刷する前に追加します（C（ "受け入れられた割合はベータ版のために描く："、100 * n.accept /イテレーション））受け入れられたとして、それはあなたを助けている場合、この答えをマーク

してください。