[PYTHON] Beispiel Adaptive MCMC

Overview Ich habe in Proseeding of NIPS2019 eine neue MCMC namens "Sample Adaptive MCMC $ ^ {[1]} $" gefunden, also habe ich es versucht. Bereits in Numpyro implementiert. Wenn Sie interessiert sind, schauen Sie sich bitte das Papier an. Die Anpassung der vorgeschlagenen Verteilung ist bei der MCMC-Methode wichtig, aber diese Methode ist eine Methode zur adaptiven Anpassung der Parameter (Mittelwert, Kovarianzmatrix) der vorgeschlagenen parametrischen Verteilung (Normalverteilung usw.) basierend auf den Stichprobenwerten. Darüber hinaus ist es eine Methode, die es einfach macht, eine effektive Stichprobengröße zu erhalten. 　[1] Michael Zhu, Sample Adaptive MCMC (https://papers.nips.cc/paper/9107-sample-adaptive-mcmc)
Algorithm Der Umriss des Algorithmus ist neu aus der vorgeschlagenen parametrischen Verteilung, die den Durchschnitt $ μ $ und die Kovarianzmatrix $ Σ $ aus den Stichprobenwerten von $ N $ berechnet und den Durchschnitt $ μ $ und die Kovarianzmatrix $ Σ $ als Parameter verwendet. Abtasten Sie den Abtastwert $ θ_ {N + 1} $. Verwerfen Sie wahrscheinlich einen dieser $ N + 1 $ -Stichprobenwerte und berechnen Sie den Durchschnitt von $ μ $ und der Kovarianzmatrix $ Σ $ der neu konstruierten $ N $ -Stichprobenwerte. Die Prozedur von ... wird wiederholt, und der Summensatz der schließlich erhaltenen Abtastwerte ist eine Annäherung an die Zielverteilung. Der Algorithmus ist wie folgt. Hier ist $ θ $, das das Probenziel ist, ein mehrdimensionaler Vektor.

__ * Eingabe * __: Objektive Verteilung $ p (θ) $, anfängliche Angebotsverteilung $ q_0 (・) $, parametrische Angebotsverteilung $ q (・ | μ, Σ) $, Anzahl der Partikel $ N $, Anzahl der Proben $ K $ , Einbrennzahl $ κ $ 　　Step.1 Holen Sie sich $ N $ Beispielwerte $ S \ _0 = $ {$ θ \ _1, θ \ _2, ... θ \ _N $} aus der ursprünglichen Angebotsverteilung $ q \ _0 (・) $. Sei $ k = 1 $. 　　Step.2

Ersetzen Sie $ S \ _ {k-1} $ durch $ S $ und berechnen Sie den Durchschnitt von $ μ $ und die Kovarianzmatrix $ Σ $ aus dem Stichprobenwert $ S $.
Erhalten Sie einen neuen Stichprobenwert $ θ \ _ {N + 1} $ aus der parametrischen Angebotsverteilung $ q (・ | μ, Σ) $ und addieren Sie ihn zu $ S $.
Erstellen Sie eine $ N + 1 $ Menge von $ S \ _ {-n} $ ohne den $ n $ -ten Stichprobenwert $ θ \ _n $ von $ S $ und teilen Sie ihn mit dem Durchschnitt $ μ \ _ {-n} $ Berechnen Sie $ N + 1 $ der Verteilungsmatrix $ Σ \ _ {-n} $. ・ $ Λ \ _n = q (θ \ _n | μ \ _ {-n}, Σ \ _ {-n}) / p (θn) $ für alle $ n (= 1 bis N + 1) $ Und standardisieren Sie so, dass die Summe $ 1 $ beträgt. ・ Als die Wahrscheinlichkeit, dass $ λ \ _n $ $ z \ _n = 1 $ ist, erhalten Sie $ n $, das $ z \ _n = 1 $ ist, aus der kategorialen Verteilung und $ S \ _k = S \ _ {-n} Lass es $ sein. 　　Step.3 Wiederholen Sie __ * Schritt 2 * __, bis die Anzahl der Wiederholungen die Anzahl der Samples $ K $ erreicht. __ * Ausgabe * __: Summenmenge des Elements $ θ $ von $ S_k (k = κ + 1 ～ K) $

Die folgende Abbildung ist ein Bilddiagramm, wenn $ N = 3 $ im Papier beschrieben ist.

Wenn $ θ $ mehrdimensional ist, ist die Berechnungslast hoch, wenn eine Kovarianzmatrix verwendet wird. Daher ist es möglich, die Geschwindigkeit zu erhöhen, indem die nicht diagonalen Komponenten der Kovarianzmatrix auf $ 0 $ angenähert werden. Es ist auch möglich, eine gemischte Gaußsche Verteilung für die vorgeschlagene parametrische Verteilung zu verwenden.

Implementation Ich habe es in Python implementiert, um das Verständnis zu fördern. Unter der Annahme eines linearen Modells werden hier auch die Planungsmatrix, die beobachteten Werte und die Messgenauigkeit als Eingaben angegeben.

# SA-MCMC
import numpy as np
from numba import jit

#Log Wahrscheinlichkeitsdichte der nicht standardisierten Zielverteilung 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

#Log-Wahrscheinlichkeitsdichte der nicht standardisierten vorgeschlagenen Verteilung q
@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

#Stichprobe aus der vorgeschlagenen Verteilung 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

#Sequentielle Aktualisierung des Durchschnittswerts
@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      :Planungsmatrix
    y      :gemessener Wert
    ram    :Meßgenauigkeit
    N      :Anzahl der Partikel
    K      :Anzahl von Beispielen
    paradim:Geschätzte Anzahl der Parameter
    U      :Einheitliche Zufallszahl(0～1),Wird aufgrund von Einschränkungen durch numba als Argument verwendet
    """
    # Step.1:Initialisieren
    #Variablendefinition
    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)
    
    #Einstellen der Anfangswerte der logarithmischen Wahrscheinlichkeitsdichte, die aus dem Mittelwert, der Kovarianzmatrix und der vorgeschlagenen Verteilung erhalten wurden
    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:Probenahme
    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)
        
        #Ablehnungsschritt:J=N+Angenommen, wenn nicht 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:Entspricht einem Summensatz
    chian = chain[:cnt-1, :]
    return chain

Experiment In diesem Experiment wird die Bayes'sche lineare Regression angestrebt. (a) Erstellen experimenteller Daten

# 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

Die Parameterschätzungsergebnisse sind wie folgt (Burnin: 2000). Die effektive Stichprobengröße (ess) ist sehr groß. Die Wahrscheinlichkeitsbewertung durch die Zielverteilung erfolgt einmal pro Stichprobe. Wenn jedoch die Anzahl der Partikel $ N $ und die Parameterdimension $ M $ beträgt, wird die Berechnung von $ O (NM ^ 2) $ zur Wahrscheinlichkeitsbewertung hinzugefügt. Es braucht Zeit. Es ist möglich, die zusätzliche Berechnungszeit auf $ O (NM) $ zu reduzieren, indem die nicht diagonale Komponente der Kovarianzmatrix auf 0 gesetzt wird. In diesem Fall scheint die Leistung jedoch erheblich verringert zu sein.

新しいビットマップイメージ.png

Summary Ich habe Sample Adaptive MCMC implementiert. Wie war es?