單組分大都市，黑斯廷斯

所以，讓我們說我有，我想從（二元正態分佈的混合物）來樣以下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

和下面的建議分佈（二元隨機遊走） -

def T(x,y,sigma): 

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

下面是在每次迭代更新「整個」狀態的都市黑斯廷斯碼 -

#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

但是，我想更新「每個組件」（即每個組件）。組件式更新）。有沒有一個簡單的方法來做到這一點？

謝謝您的幫助！

來源

2017-07-18 tattybojangler

您首先必須計算目標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））''） – misterkugelblitz

從你的問題的語氣中，我假設你正在尋找性能改進。 MonteCarlo算法計算密集。如果你在比python這樣的解釋型語言更低級別的算法中執行，你會得到更好的結果。寫一個c-extension。

也有實現可用於Python（PyStan，PyMC3）。

來源

2017-07-18 22:17:59

單組分大都市，黑斯廷斯

回答

相關問題