单组件 Metropolis-Hastings答案

【问题标题】：Single Component Metropolis-Hastings单组件 Metropolis-Hastings
【发布时间】：2017-07-18 20:43:49
【问题描述】：

所以，假设我有以下二维目标分布，我想从中采样（二元正态分布的混合） -

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

以及以下提案分布（二元随机游走）-

def T(x,y,sigma):

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

以下是 Metropolis Hastings 代码，用于在每次迭代中更新“整个”状态 -

#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

但是，我想在每次迭代中更新“每个组件”（即逐个组件更新）。有没有一种简单的方法可以做到这一点？

感谢您的帮助！

【问题讨论】：

您首先必须计算目标 PDF 的边际 PDF。然后您可以对组件进行采样 Y[i] = X[i]+stats.norm.rvs(0,sigma,1) 并接受/拒绝组件（即 r = (marg_targ_dist(Y[i])*T(Y[i],X[i],sigma))/(marg_targ_dist(X[i])*T(X[i],Y[i],sigma)) ）

标签： python statistics montecarlo markov-chains mcmc

【解决方案1】：

从您的问题的语气来看，我认为您正在寻找性能改进。蒙特卡洛算法是计算密集型的。如果您在比 Python 等解释性语言更低级别的算法中执行，您将获得更好的结果，例如写一个c扩展。

还有可用于 python 的实现（PyStan、PyMC3）。

【讨论】：