PRML §3.3.1 Reproduire le diagramme de convergence de la distribution des paramètres par régression linéaire bayésienne

Si vous disposez de PRML et de planches couleur qui sont promus lors de sessions d'étude internes, vous voudrez les reproduire sérieusement. Cette fois, j'ai reproduit l'illustration de la convergence de la probabilité postérieure de la distribution des paramètres en 3.3.1. Ce chiffre est

  1. Paramètre de vraisemblance basé sur les données d'observation les plus récentes
  2. Distribution de probabilité pré / post des paramètres
  3. Ajustement des résultats avec les paramètres obtenus à partir de la distribution de probabilité pré / post

Il est intéressant de voir comment converge la probabilité ex post facto.

supposition

Pièces de calcul requises

Création d'une matrice de planification $ \ Phi $ à partir des données d'observation

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

Fonction pour trouver la moyenne $ m_N = \ beta S_N \ Phi ^ {T} t $ après avoir observé les données N fois

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)

Fonction pour trouver la covariance $ S_N = (\ alpha I + \ beta \ Phi ^ {T} \ Phi) ^ {-1} $ après avoir observé les données N fois

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))

Espérance des données d'observation

Puisque les données observées ont un bruit gaussien avec une précision $ \ beta $, il s'agit d'une fonction de vraisemblance logarithmique.

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

Utiliser. Trouvez $ w $ dans la plage que vous souhaitez tracer avec la fonction de vraisemblance du journal.

def calc_likelifood(beta, t, x, w):
    """
Trouvez la vraisemblance logarithmique d'une observation
    """
    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):
    """
Graphique de la probabilité de la valeur observée
    """
    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)

Probabilité post-hoc

La distribution normale bivariée avec $ m_N $ et $ S_N $ comme paramètres obtenus par la formule ci-dessus est la distribution de probabilité postérieure de $ w $. Dans le cas de $ m_0 $ et $ S_0 $, la distribution de probabilité antérieure.

def plot_probability(mean, cov, title='', ax=None):
    """
Distribution de probabilité(Gauss bivarié)Terrain
    """
    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)

Espace de données

Échantillonnez $ w $ 6 fois à partir de la distribution normale bivariée avec $ m_N $ et $ S_N $ comme paramètres obtenus par la formule ci-dessus et tracez-la.

résultat

Avant l'échantillon de données index.png

Après un échantillon de données index2.png

Après 8 échantillons de données index3.png

C'est fait. Tout le code est également posté sur github. https://github.com/hagino3000/public-ipynb/blob/master/PRML/PRML%203.3.ipynb

Recommended Posts

PRML §3.3.1 Reproduire le diagramme de convergence de la distribution des paramètres par régression linéaire bayésienne
Reproduction sur plaque de régression linéaire bayésienne (PRML §3.3)
Implémentation python de la classe de régression linéaire bayésienne
Implémenté en Python PRML Chapitre 3 Régression linéaire bayésienne
Diagramme PRML Figure 1.15 Biais dans l'estimation la plus probable de la distribution gaussienne