monocomposant Metropolis-Hastings

Question

Alors, disons que je la répartition suivante cible 2 dimensions que je voudrais échantillonner (un mélange de distributions bivariées normales) -

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 

et la distribution de proposition suivante (une marche aléatoire bivariée) -

def T(x,y,sigma): 

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

Voici le code Metropolis Hastings pour mettre à jour l'état "entier" dans chaque itération -

#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 

Cependant, je voudrais mettre à jour « par composant » (à savoir mise à jour par composant) à chaque itération. Y a-t-il un moyen simple de le faire?

Nous vous remercions de votre aide!

Answer 1

Du ton de votre question, je suppose que vous êtes à la recherche des améliorations de performance. Les algorithmes de MonteCarlo demandent beaucoup de ressources. Vous obtiendrez de meilleurs résultats si vous effectuez des algorithmes à un niveau inférieur à celui d'un langage interprété comme python, par ex. écrire une c-extension.

Il existe également des implémentations disponibles pour python (PyStan, PyMC3).