[PYTHON] PRML §3.3.1 Reproduzieren Sie das Konvergenzdiagramm der Parameterverteilung durch Bayes'sche lineare Regression

Wenn Sie PRML- und Farbtafeln haben, die bei internen Studiensitzungen beworben werden, sollten Sie diese ernsthaft reproduzieren. Dieses Mal habe ich die Illustration reproduziert, wie die hintere Wahrscheinlichkeit der Verteilung der Parameter in 3.3.1 konvergiert. Diese Zahl ist

1. Parameterwahrscheinlichkeit basierend auf den neuesten Beobachtungsdaten
2. Verteilung der Parameter vor / nach der Wahrscheinlichkeit
3. Anpassung der Ergebnisse mit Parametern aus der Wahrscheinlichkeitsverteilung vor / nach

Es ist interessant zu sehen, wie die Ex-post-Wahrscheinlichkeit konvergiert.

Annahme

--Betrachten Sie eine eindimensionale Eingabevariable $ x $ und eine eindimensionale Zielvariable $ t $ -Fit zu einem linearen Modell von $ y (x, w) = w_0 + w_1x $ --Modellparameter Die beiden Verteilungen sind bivariate Gaußsche Verteilungen.

• Rauschen ist Gaußsches Rauschen

Erforderliche Berechnungsteile

Erstellen einer Planungsmatrix $ \ Phi $ aus Beobachtungsdaten

def design_matrix(x):
    return np.array([[1, xi] for xi in x])

Funktion zum Ermitteln des Durchschnitts $ m_N = \ beta S_N \ Phi ^ {T} t $ nach N-maliger Beobachtung der Daten

def calc_mn(alpha, beta, x, t):
    Phi = design_matrix(x)
    Sn = calc_Sn(alpha, beta, x)
    return beta * Sn.dot(Phi.T).dot(t)

Funktion zum Ermitteln der Kovarianz $ S_N = (\ alpha I + \ beta \ Phi ^ {T} \ Phi) ^ {-1} $ nach N-maliger Beobachtung der Daten

I = np.identity(2)

def calc_Sn(alpha, beta, x):
    Phi = design_matrix(x)
    return np.linalg.inv(alpha*I + beta*Phi.T.dot(Phi))

Erwartbarkeit der Beobachtungsdaten

Da die beobachteten Daten ein Gaußsches Rauschen mit der Genauigkeit $ \ beta $ aufweisen, handelt es sich um eine logarithmische Wahrscheinlichkeitsfunktion.

logL(w) = -\frac{\beta}{2}(t-w^{T}\phi(x))^2 + cons

Benutzen. Suchen Sie $ w $ in dem Bereich, den Sie mit der Protokollwahrscheinlichkeitsfunktion zeichnen möchten.

def calc_likelifood(beta, t, x, w):
    """
Finden Sie die logarithmische Wahrscheinlichkeit einer Beobachtung
    """
    w = np.array(w)
    phi_x = np.array([1, x])
    return -1 * beta / 2 * (t - w.T.dot(phi_x))**2

def plot_likelifood(beta, t, x, title='', ax=None):
    """
Darstellung der beobachteten Wertwahrscheinlichkeit
    """
    w0 = np.linspace(-1, 1, 100)
    w1 = np.linspace(-1, 1, 100)
    W0,W1 = np.meshgrid(w0, w1)
    L = []
    for w0i in w0:
        L.append([calc_likelifood(beta, t, x, [w0i, w1i]) for w1i in w1])

    ax.pcolor(W0, W1, np.array(L).T, cmap=plt.cm.jet, vmax=0, vmin=-1)
    ax.set_xlabel('$w_0$')
    ax.set_ylabel('$w_1$')
    ax.set_title(title)

Post-hoc-Wahrscheinlichkeit

Die bivariate Normalverteilung mit $ m_N $ und $ S_N $ als Parametern, die durch die obige Formel erhalten werden, ist die hintere Wahrscheinlichkeitsverteilung von $ w $. Im Fall von $ m_0 $ und $ S_0 $ die Wahrscheinlichkeitsverteilung.

def plot_probability(mean, cov, title='', ax=None):
    """
Wahrscheinlichkeitsverteilung(Bivariate Gauß)Handlung
    """
    w0 = np.linspace(-1, 1, 100)
    w1 = np.linspace(-1, 1, 100)
    W0,W1 = np.meshgrid(w0, w1)
    P = []
    for w0i in w0:
        P.append([scipy.stats.multivariate_normal.pdf([w0i,w1i], mean, cov) for w1i in w1])
        
    ax.pcolor(W0, W1, np.array(P).T, cmap=plt.cm.jet)
    ax.set_xlabel('$w_0$')
    ax.set_ylabel('$w_1$')
    ax.set_title(title)

Datenraum

Probieren Sie $ w $ 6-mal aus der bivariaten Normalverteilung mit $ m_N $ und $ S_N $ als Parametern, die durch die obige Formel erhalten wurden, und zeichnen Sie sie auf.

Ergebnis

Vor der Datenprobe

Nach einer Datenprobe

Nach 8 Datenproben

Es ist fertig. Der gesamte Code wird auch auf github veröffentlicht. https://github.com/hagino3000/public-ipynb/blob/master/PRML/PRML%203.3.ipynb