Implémenté en Python PRML Chapitre 1 Estimation bayésienne

Dans la continuité de la dernière fois, il s'agit d'une implémentation de l '"estimation bayésienne 1.2.6" du chapitre 1 de PRML. Il s'agit d'une reproduction de la figure 1.17, mais avec M = 9, la moyenne de la distribution prédite et la plage de ± 1σ sont indiquées.

C'était abstrait quand j'ai suivi la formule, et cela ne m'est pas venu. En l'implémentant, j'ai été frappé par la visualisation concrète que (1.70) et (1.71) montrent la distribution. Il ne peut pas être déterminé à partir de cet exemple en quoi cette précision de prédiction diffère de l'ajustement de courbe polymorphe précédent, Cela a été très utile pour comprendre à quoi chacun ressemble.

Flux de mise en œuvre approximatif

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

 p(t| x, {\bf x}, {\bf t}) = N(t| m(x), s^2(x)) (1.69)

② Cette fois, la fonction de base est $ \ phi_i (x) = x ^ i $ pour $ i = 0, ... M $.

(3) Puisqu'il existe les trois formules suivantes pour calculer la moyenne et la variance de la distribution prédite, implémentez chacune d'elles pour visualiser la distribution dans la plage prédite de 0,0 $ <x <1,0 $.

m(x) = \beta\phi(x)^{\bf T}{\bf S} \sum_{n=1}^N \phi(x_n)t_n(1.70)]

 s^2(x) = \beta^{-1} + \phi(x)^{\bf T} {\bf S} \phi(x)(1.71)

{\bf S}^{-1} = \alpha {\bf I} + \beta \sum_{n=1}^N \phi(x_n)\phi(x_n)^{\bf T}(1.72)

code

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

#From p31 the auther define phi as following
def phi(x):
    return np.array([x**i for i in xrange(M+1)]).reshape((M+1, 1))

#(1.70) Mean of predictive distribution
def m(x, x_train, y_train, S):
    sum = np.array(zeros((M+1, 1)))
    for n in xrange(len(x_train)):
        sum += np.dot(phi(x_train[n]), y_train[n])
    return Beta * phi(x).T.dot(S).dot(sum)

    
#(1.71) Variance of predictive distribution   
def s2(x, S):
    return 1.0/Beta + phi(x).T.dot(S).dot(phi(x))

#(1.72)
def S(x_train, y_train):
    I = np.identity(M+1)
    Sigma = np.zeros((M+1, M+1))
    for n in range(len(x_train)):
        Sigma += np.dot(phi(x_train[n]), phi(x_train[n]).T)
    S_inv = alpha*I + Beta*Sigma
    S = inv(S_inv)
    return S

if __name__ == "__main__":
    alpha = 0.005
    Beta = 11.1
    M = 9
    
    #Sine curve
    x_real = np.arange(0, 1, 0.01)
    y_real = np.sin(2*np.pi*x_real)
    
    ##Training Data
    N=10
    x_train = np.linspace(0, 1, 10)
    
    #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,N)
    
    
    S = S(x_train, y_train)
    
    #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 1.17")

résultat