RのMetropolis-Hastingsアルゴリズム：正しい結果はありますか？

Question

私のメトロポリス - ヘイスティングス問題は固定の二項分布を持ち、すべての提案分布q（i、j）は0.5です。プロットとヒストグラムを参照して、アルゴリズムを0.5付近に集中させると、二項分布からの確率は？

pi <- function(x, n, p){ 
# returning value from binomial distribution 
    result <- (factorial(n)/(factorial(x) * factorial(n - x))) * 
     p^x * (1 - p)^(n - x) 
    return(result) 
} 


metropolisAlgorithm <- function(n, p, T){ 
# implementation of the algorithm 
# @n,p binomial parameters 
# @T number of time steps 
    X <- rep(runif(1),T) 
    for (t in 2:T) { 
     Y <- runif(1) 
     alpha <- pi(X[t - 1], n, p)/pi(Y, n, p) 
     if (runif(1) < alpha) X[t] <- Y 
     else X[t] < X[t - 1] 
    } 
    return(X) 
} 

# calling M-H algorithm and plotting result 
test <- metropolisAlgorithm(40,0.5,5000) 
par(mfrow=c(2,1)) 
plot(test, type = "l") 
hist(test, breaks = 40) 

Answer 1

あなたは3つの問題だった：あなたのランダムウォークが範囲1:nの整数ではなく実数でかけなければなりませんので、あなたは、二項分布をシミュレートするように見える

1）範囲は[0,1]です。

2）あなたは分子を持っていたし、分母はあなたがX[t] < X[t - 1]にタイプミスがあったalpha

3）の計算に切り替えます。

（@ZheyuanLiが提案されているようdbinom機能の使用を含む）ビットをこれらを修正し、あなたのコードをクリーンアップ利回り：

metropolisAlgorithm <- function(n, p, T){ 
    # implementation of the algorithm 
    # @n,p binomial parameters 
    # @T number of time steps 
    X <- rep(0,T) 
    X[1] <- sample(1:n,1) 
    for (t in 2:T) { 
    Y <- sample(1:n,1) 
    alpha <- dbinom(Y,n,p)/dbinom(X[t-1],n,p) 
    if (runif(1) < alpha){ 
     X[t] <- Y 
    }else{ 
     X[t] <- X[t - 1]} 
    } 
    return(X) 
} 

# calling M-H algorithm and plotting result 
test <- metropolisAlgorithm(40,0.5,5000) 
par(mfrow=c(2,1)) 
plot(test, type = "l") 
hist(test) breaks = 40) 

（完璧な理にかなっている）の典型的な出力：