[PYTHON] 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

Paramètre de vraisemblance basé sur les données d'observation les plus récentes
Distribution de probabilité pré / post des paramètres
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

Considérez une variable d'entrée unidimensionnelle $ x $ et une variable cible unidimensionnelle $ t $ -Fit à un modèle linéaire de $ y (x, w) = w_0 + w_1x $ --Paramètres du modèle Les deux distributions sont des distributions gaussiennes bivariées. --Le bruit est un bruit gaussien

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

Après un échantillon de données

Après 8 échantillons de données

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