Implémenté en Python PRML Chapitre 3 Régression linéaire bayésienne

Reconnaissance de formes et apprentissage automatique, le chapitre 3 reproduit comment la régression linéaire bayésienne fait converger la distribution prédite. Ce que nous faisons est presque la même que l'estimation bayésienne que nous avons faite au chapitre 1, seule la fonction de base est la fonction gaussienne (3.4). Par conséquent, le code est presque détourné.

Par rapport au chapitre 1, au chapitre 3, l'expression de la formule mathématique est réécrite pour manipuler les vecteurs, y a-t-il une intention dans ce domaine? Depuis que j'ai introduit le concept de vecteur de caractéristiques dans (3.16) ici, je me demande si c'est la manière de l'écrire, mais qu'en est-il?

Flux de mise en œuvre approximatif

(1) Ce que vous voulez trouver est la distribution prévue (3,57).

 p(t|{\bf x}, {\bf t}, \alpha, \beta) = N (t| {\bf m}^T_N \phi({\bf x}),  \sigma^2_N ({\bf x})) (3.58)

② La fonction de base est $ \ phi_i (x) = \ exp \ {- \ frac {(x- \ mu_j) ^ 2} {2s ^ 2} } $ in (3.4). Bien que 9 fonctions de base soient spécifiées, elles sont organisées de manière à être uniformément réparties car il n'y a aucune instruction sur la façon de les organiser.

(3) Puisqu'il existe les trois formules suivantes pour calculer la moyenne et la variance de la distribution prédite, appliquez chacune d'elles pour obtenir la distribution prédite.

{\bf m}_N = \beta{\bf S}_N\Phi(x_n)^{\bf T}{\bf t}(3.53)

\sigma^2_N(x) = \beta^{-1} + \phi({\bf x})^{\bf T} {\bf S}_N \phi({\bf x}). (3.59)

{\bf S}^{-1}_N = \alpha {\bf I} + \beta\Phi^{\bf T}\Phi(3.54)

code

En dessinant la figure, @ naoya_t's Cet article (Reproduction of the plate of bayesian linear regression (PRML§3.3)) montre comment les données de l'échantillon convergent à mesure qu'elles augmentent. ) Je n'étais pas convaincu que je pouvais trouver une idée pour dessiner plus joliment, alors je l'ai utilisée comme référence, presque comme une copie. Merci beaucoup.

import numpy as np
from numpy.linalg import inv
import pandas as pd
from pylab import *
import matplotlib.pyplot as plt

def func(x_train, y_train):
    # (3.4) Gaussian basis function
    def phi(s, mu_i, M, x):
        return np.array([exp(-(x - mu)**2 / (2 * s**2)) for mu in mu_i]).reshape((M, 1))

    #(3.53)' ((1.70)) Mean of predictive distribution
    def m(x, x_train, y_train, S):
        sum = np.array(zeros((M, 1)))
        for n in xrange(len(x_train)):
            sum += np.dot(phi(s, mu_i, M, x_train[n]), y_train[n])
        return Beta * phi(s, mu_i, M, x).T.dot(S).dot(sum)
    
    #(3.59)'((1.71)) Variance of predictive distribution   
    def s2(x, S):
        return 1.0/Beta + phi(s, mu_i, M, x).T.dot(S).dot(phi(s, mu_i, M, x))

    #(3.53)' ((1.72))
    def S(x_train, y_train):
        I = np.identity(M)
        Sigma = np.zeros((M, M))
        for n in range(len(x_train)):
            Sigma += np.dot(phi(s, mu_i, M, x_train[n]), phi(s, mu_i, M, x_train[n]).T)
        S_inv = alpha*I + Beta*Sigma
        S = inv(S_inv)
        return S
    
    #params for prior probability
    alpha = 0.1
    Beta = 9
    s = 0.1

    #use 9 gaussian basis functions
    M = 9

    # Assign basis functions
    mu_i = np.linspace(0, 1, M)
    
    S = S(x_train, y_train)

    #Sine curve
    x_real = np.arange(0, 1, 0.01)
    y_real = np.sin(2*np.pi*x_real)
    
    #Seek predictive distribution corespponding to entire x
    mean = [m(x, x_train, y_train, S)[0,0] for x in x_real]
    variance = [s2(x, S)[0,0] for x in x_real]
    SD = np.sqrt(variance)
    upper = mean + SD
    lower = mean - SD
    
    plot(x_train, y_train, 'bo')
    plot(x_real, y_real, 'g-')
    plot(x_real, mean, 'r-')
    fill_between(x_real, upper, lower, color='pink')
    xlim(0.0, 1.0)
    ylim(-2, 2)
    title("Figure 3.8")
    show()

if __name__ == "__main__":
    # Maximum observed data points (underlining function plus some noise)
    x_train = np.linspace(0, 1, 26)

    #Set "small level of random noise having a Gaussian distribution"
    loc = 0
    scale = 0.3
    y_train =  np.sin(2*np.pi*x_train) + np.random.normal(loc,scale,26)


    #Sample data pick up
    def randidx(n, k):
        r = range(n)
        shuffle(r)
        return sort(r[0:k])

    for k in (1, 2, 5, 26):
        indices = randidx(size(x_train), k)
        func(x_train[indices], y_train[indices])

résultat

référence

Reproduction de la régression linéaire bayésienne (PRML §3.3)