[PYTHON] Exemple de MCMC adaptatif
- Overview J'ai trouvé un nouveau MCMC appelé "Sample Adaptive MCMC $ ^ {[1]} $" dans Proseeding de NIPS2019, alors je l'ai essayé. Déjà implémenté dans Numpyro. Si vous êtes intéressé, jetez un œil à l'article. L'ajustement de la distribution proposée est important dans la méthode MCMC, mais cette méthode est une méthode d'ajustement adaptatif des paramètres (moyenne, matrice de covariance) de la distribution paramétrique proposée (distribution normale, etc.) en fonction des valeurs de l'échantillon. De plus, c'est une méthode qui permet d'obtenir facilement une taille d'échantillon efficace. [1] Michael Zhu, Sample Adaptive MCMC (https://papers.nips.cc/paper/9107-sample-adaptive-mcmc)
- Algorithm Le contour de l'algorithme est nouveau à partir de la distribution paramétrique proposée qui calcule la moyenne $ μ $ et la matrice de covariance $ Σ $ à partir des valeurs d'échantillon de $ N $ et utilise la moyenne $ μ $ et la matrice de covariance $ Σ $ comme paramètres. Échantillon de la valeur d'échantillon $ θ_ {N + 1} $. Supprimez probablement l'une de ces valeurs d'échantillon $ N + 1 $ et calculez la moyenne $ μ $ et la matrice de covariance $ Σ $ des valeurs d'échantillon $ N $ nouvellement construites. La procédure de ... est répétée, et l'ensemble de somme des valeurs d'échantillon finalement obtenues est une approximation de la distribution cible. L'algorithme est le suivant. Ici, $ θ $, qui est la cible de l'échantillon, est un vecteur multidimensionnel.
__ * Entrée * __: Distribution objective $ p (θ) $, distribution de proposition initiale $ q_0 (・) $, distribution de proposition paramétrique $ q (・ | μ, Σ) $, nombre de particules $ N $, nombre d'échantillons $ K $ , Numéro de rodage $ κ $ Step.1 Obtenez des exemples de valeurs $ N $ $ S \ _0 = $ {$ θ \ _1, θ \ _2, ... θ \ _N $} à partir de la distribution de proposition initiale $ q \ _0 (・) $. Soit $ k = 1 $. Step.2 -Substituer $ S \ _ {k-1} $ dans $ S $ et calculer la moyenne $ μ $ et la matrice de covariance $ Σ $ à partir de la valeur d'échantillon $ S $. -Obtenir une nouvelle valeur d'échantillon $ θ \ _ {N + 1} $ à partir de la distribution de proposition paramétrique $ q (・ | μ, Σ) $ et l'ajouter à $ S $. -Créez un ensemble $ N + 1 $ de $ S \ _ {-n} $ excluant la valeur d'échantillon $ n $ ème $ θ \ _n $ de $ S $, et partagez-le avec la moyenne $ μ \ _ {-n} $ Calculez $ N + 1 $ de la matrice de distribution $ Σ \ _ {-n} $. ・ $ Λ \ _n = q (θ \ _n | μ \ _ {-n}, Σ \ _ {-n}) / p (θn) $ pour tout $ n (= 1 à N + 1) $ Et standardisez pour que la somme soit de 1 $. -Obtient $ n $ de la distribution catégorielle comme la probabilité que $ λ \ _n $ soit $ z \ _n = 1 $, et $ S \ _k = S \ _ {-n} Que ce soit $. Step.3 Répétez __ * Étape 2 * __ jusqu'à ce que le nombre de répétitions atteigne le nombre d'échantillons $ K $. __ * Sortie * __: Somme des éléments $ θ $ de $ S_k (k = κ + 1 ～ K) $
La figure ci-dessous est un diagramme d'image lorsque $ N = 3 $ décrit dans l'article.
Si $ θ $ est multidimensionnel, la charge de calcul sera élevée si une matrice de covariance est utilisée. Par conséquent, il est possible d'augmenter la vitesse en approchant les composantes hors diagonale de la matrice de covariance à $ 0 $. Il est également possible d'utiliser une distribution gaussienne mixte pour la distribution paramétrique proposée.
- Implementation Je l'ai implémenté en Python pour favoriser la compréhension. Ici, en supposant un modèle linéaire, la matrice de planification, les valeurs observées et la précision des mesures sont également données comme entrées.
# SA-MCMC
import numpy as np
from numba import jit
#Densité de probabilité logarithmique de la distribution cible non normalisée p
@jit("f8(f8[:,:],f8[:],f8[:],f8)", nopython=True)
def target(X, y, w, ram):
mu_w0 = 0.0
ram_w0 =1.0
logp1 = np.sum(-0.5 * ram_w0 * (w - mu_w0) ** 2)# logP(w)+C
logp2 = np.sum(-0.5 * ram * (y - X @ w) ** 2) # logP(y|X, w, ram)+C
logp = logp1 + logp2
return logp
#Densité de probabilité logarithmique de la distribution proposée non standardisée
@jit("f8(f8[:],f8[:,:],f8[:])", nopython=True)
def proposal(mu, cov, sita):
A = -0.5 * np.linalg.slogdet(cov)[1]
B = -0.5 * (sita - mu).T @ np.linalg.solve(cov, (sita - mu))
return A+B
#Échantillonnage à partir de la distribution proposée q
@jit("f8[:](f8[:],f8[:,:],i8)",nopython=True)
def sampling(mu, cov, w_dim):
L = np.linalg.cholesky(cov)
t = np.random.randn(w_dim)
return L @ t + mu
#Mise à jour séquentielle de la valeur moyenne
@jit("f8[:](f8[:],f8[:],f8[:],i8)",nopython=True)
def update_mean(mean, new_s, old_s, N):
mu = mean + 1/N * (new_s - old_s)
return mu
@jit("f8[:,:](f8[:,:],f8[:],f8,i8,i8,i8,f8[:])",nopython=True)
def SA_MCMC(X, y, ram, N, K, paradim, U):
"""
X :Matrice de planification
y :la valeur de mesure
ram :Précision de mesure
N :Nombre de particules
K :Le nombre d'échantillons
paradim:Nombre estimé de paramètres
U :Nombre aléatoire uniforme(0～1),Utilisé comme argument en raison des restrictions de numba
"""
# Step.1:Initialisation
#Définition variable
np.random.seed(8)
chain = np.zeros((K, paradim))
S = np.random.randn(N+1, paradim)
S_mean = np.zeros(paradim)
S_replace = np.zeros((N+1, paradim))
q = np.zeros(N+1)
r = np.zeros(N+1)
#Définition des valeurs initiales de la densité de probabilité logarithmique obtenue à partir de la moyenne, de la matrice de covariance et de la distribution proposée
for d in range(paradim):
S_mean[d] = S[:-1, d].mean()
S_cov = np.cov(S[:-1, :].T)
for n in range(N):
q[n] = target(X, y, S[n, :], ram)
# Step.2:échantillonnage
cnt = 0
for i in range(K):
s = sampling(S_mean, S_cov, paradim)
q[-1] = target(X, y, s, ram)
for n in range(N+1):
S_replace[:, :] = S[:, :]
S_replace[n, :] = s
S_mean_replace = update_mean(S_mean, s, S[n,:], N)
S_cov_replace = np.cov(S_replace[:-1, :].T)
p = proposal(S_mean_replace, S_cov_replace, sita=S[n, :])
r[n] = p - q[n]
rmax = r.max()
r = np.exp(r - rmax)
r = r / np.sum(r)
#Étape de rejet:J=N+Adopté si différent de 1
c = 0
u=U[i]
for j in range(N+1):
c += r[j]
if c > u:
J = j
break
S_mean = update_mean(S_mean, s, S[J,:], N)
S[J, :] = s
S_cov = np.cov(S[:-1, :].T)
q[J] = q[-1]
if J != N+1:
chain[cnt] = s
cnt += 1
# Step.3:Équivalent à un ensemble de somme
chian = chain[:cnt-1, :]
return chain
- Experiment Dans cette expérience, la régression linéaire bayésienne est ciblée. (a) Création de données expérimentales
# make data
def makedata(n, w_dim, vmin, vmax, seed, ram=0.0, Train=False):
np.random.seed(seed)
x = np.linspace(vmin, vmax, n)
w = np.random.uniform(-1, 1, w_dim)
X = np.zeros((n, w_dim))
for d in range(w_dim):
X[:, d] = x ** d
y = X @ w
if Train:
y += np.random.randn(n) / ram ** 0.5
return x, w, X, y
w_dim = 4
vmin = -2
vmax = 2
seed = 29
ram = 5
x_true, w_true, X_true, y_true = makedata(2000, w_dim, vmin-2, vmax+2, seed, Train=False)
x_train, _, X_train, y_train = makedata(30, w_dim, vmin, vmax, seed, ram, Train=True)
print(w_true)
# [ 0.72751997 -0.43018807 -0.85348722 0.52647441]
(b) SA-MCMC
%%time
sample_n = 10000
U1 = np.random.uniform(size=sample_n)
chain1 = SA_MCMC(X_train, y_train, ram, 100, sample_n, w_dim, U1)
U2 = np.random.uniform(size=sample_n)
chain2 = SA_MCMC(X_train, y_train, ram, 100, sample_n, w_dim, U2)
U3 = np.random.uniform(size=sample_n)
chain3 = SA_MCMC(X_train, y_train, ram, 100, sample_n, w_dim, U3)
# Wall time: 19.2 s
Les résultats de l'estimation des paramètres sont les suivants (burnin: 2000). La taille effective de l'échantillon (ess) est très grande. L'évaluation de la probabilité par la distribution cible est une fois par échantillonnage, mais si le nombre de particules est $ N $ et la dimension du paramètre est $ M $, le calcul de $ O (NM ^ 2) $ est ajouté à l'évaluation de la vraisemblance. Ça prend du temps. Il est possible de réduire le temps de calcul supplémentaire à $ O (NM) $ en fixant la composante hors diagonale de la matrice de covariance à 0, mais dans ce cas, les performances semblent être considérablement réduites.
- Summary J'ai implémenté Sample Adaptive MCMC. Comment était-ce?