Implemented in Python PRML Chapter 1 Bayesian Inference

Continuing from the last time, it is an implementation of "1.2.6 Bayesian inference" from PRML Chapter 1. This is a reproduction of Figure 1.17, but with M = 9, the mean of the predicted distribution and the range of ± 1σ are shown.

It was abstract when I just followed the formula, and it didn't come to me. By implementing it, I was struck by the concrete visualization that (1.70) and (1.71) show distribution. It cannot be determined from this example how this prediction accuracy differs from the previous polynomial curve fitting, It was very helpful in understanding what each one looks like.

Rough flow of implementation

(1) What we want to find is the predicted distribution (1.69).

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

② This time, the basis function is $ \ phi_i (x) = x ^ i $ for $ i = 0, ... M $.

(3) Since there are the following three formulas for calculating the mean and variance of the predicted distribution, implement each of them to visualize the distribution in the predicted range of $ 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")

result

Screen Shot 2015-09-18 at 07.22.52.png

Recommended Posts

Implemented in Python PRML Chapter 1 Bayesian Inference
Implemented in Python PRML Chapter 3 Bayesian Linear Regression
Implemented in Python PRML Chapter 7 Nonlinear SVM
Implemented in Python PRML Chapter 1 Polynomial Curve Fitting
Implemented in Python PRML Chapter 4 Classification by Perceptron Algorithm
PRML Chapter 4 Bayesian Logistic Regression Python Implementation
PRML Chapter 1 Bayesian Curve Fitting Python Implementation
Implemented SimRank in Python
Implemented Shiritori in Python
I implemented Robinson's Bayesian Spam Filter in python
PRML Chapter 12 Bayesian Principal Component Analysis Python Implementation
A addictive point in "Bayesian inference experience with Python"
Sudoku solver implemented in Python 3
[Python] Bayesian inference with Pyro
6 Ball puzzle implemented in python
Implemented image segmentation in python (Union-Find)
100 Language Processing Knock Chapter 1 in Python
PRML Chapter 5 Neural Network Python Implementation
Bayesian optimization package GPyOpt in Python
Widrow-Hoff learning rules implemented in Python
Implemented label propagation method in Python
PRML Chapter 3 Evidence Approximation Python Implementation
Implemented Perceptron learning rules in Python
Implemented in 1 minute! LINE Notify in Python
PRML Chapter 8 Product Sum Algorithm Python Implementation
Spiral book in Python! Python with a spiral book! (Chapter 14 ~)
PRML Chapter 5 Mixed Density Network Python Implementation
A simple HTTP client implemented in Python
PRML Chapter 9 Mixed Gaussian Distribution Python Implementation
PRML Chapter 14 Conditional Mixed Model Python Implementation
PRML Chapter 10 Variational Gaussian Distribution Python Implementation
PRML Chapter 6 Gaussian Process Regression Python Implementation
PRML Chapter 2 Student's t Distribution Python Implementation
I tried using Bayesian Optimization in Python
I implemented Cousera's logistic regression in Python
Implemented Stooge sort in Python3 (Bubble sort & Quicksort)
Introduction to Effectiveness Verification Chapter 1 in Python
Duck book implemented in Python "Bayesian statistical modeling with Stan and R"
Quadtree in Python --2
Python in optimization
CURL in python
Metaprogramming in Python
Python 3.3 in Anaconda
SendKeys in Python
Introduction to effectiveness verification Chapter 3 written in Python
Epoch in Python
Discord in Python
Sudoku in Python
DCI in Python
quicksort in python
nCr in python
[Python] Implemented automation in excel file copying work
N-Gram in Python
PRML Chapter 11 Markov Chain Monte Carlo Python Implementation
Programming in python
Constant in python
Lifegame in Python.
FizzBuzz in Python
Sqlite in python
StepAIC in Python
I implemented the inverse gamma function in python