単一コンポーネントメトロポリス - ヘイスティングス

Question

それでは、私は（二変量正規分布の混合物）から採取したいと思い、次の2次元ターゲット分布を持っているとしましょう -

import numba 
import numpy as np 
import scipy.stats as stats 
import seaborn as sns 
import pandas as pd 
import matplotlib.mlab as mlab 
import matplotlib.pyplot as plt 
%matplotlib inline 

def targ_dist(x): 

target = (stats.multivariate_normal.pdf(x,[0,0],[[1,0],[0,1]])+stats.multivariate_normal.pdf(x,[-6,-6],[[1,0.9],[0.9,1]])+stats.multivariate_normal.pdf(x,[4,4],[[1,-0.9],[-0.9,1]]))/3 
return target 

と、次の提案分布（二変量ランダムウォーク） -

#Initialising 

n_iter = 30000 

# tuning parameter i.e. variance of proposal distribution 
sigma = 2 

# initial state 
X = stats.uniform.rvs(loc=-5, scale=10, size=2, random_state=None) 

# count number of acceptances 
accept = 0 

# store the samples 
MHsamples = np.zeros((n_iter,2)) 

# MH sampler 
for t in range(n_iter): 

    # proposals 
    Y = X+stats.norm.rvs(0,sigma,2) 

    # accept or reject 
    u = stats.uniform.rvs(loc=0, scale=1, size=1) 

    # acceptance probability 
    r = (targ_dist(Y)*T(Y,X,sigma))/(targ_dist(X)*T(X,Y,sigma)) 
    if u < r: 
     X = Y 
     accept += 1 
    MHsamples[t] = X 

-

def T(x,y,sigma): 

return stats.multivariate_normal.pdf(y,x,[[sigma**2,0],[0,sigma**2]]) 

次は、すべての反復で「全体」の状態を更新するためのメトロポリスヘイスティングスコードです

ただし、「コンポーネントごと」（つまり、コンポーネントごとの更新）を繰り返します。これを行う簡単な方法はありますか？

ありがとうございました！

Answer 1

あなたの質問の色調から、私はあなたがパフォーマンスの改善を見ていると思います。 MonteCarloアルゴリズムはかなり計算が集中しています。 Pythonなどのインタプリタ言語よりも低いレベルでアルゴリズムを実行すると、より良い結果が得られます。 c-extensionを書く。

Python（PyStan、PyMC3）で使用できる実装もあります。