prml_fig3_8.py

# coding:utf-8
import numpy as np
import math
import random
import matplotlib.pyplot as plt

mu = np.arange(0, 1.01, 0.125) # the positions of the Gaussian basis function
S = 0.1 # the spatial factor of the Gaussian basis function
ALPHA = 0.1 # the presition parameter of the prior
BETA = 9 # the presition parameter of the noise

def sinpi(x):
    return math.sin(2*math.pi * x)

def base(x, n):
    return np.exp(- (x - mu[n]) ** 2 / (2 * S**2))

def makeData(N):
    l = []
    for i in range(N):
        x = random.random()
        error = np.random.normal(0, 0.3)
        l.append((x, sinpi(x) + error))
    return l

def main():
    N = int(input("number of samples: "))
    M = 9
    data = makeData(N)
    phi = [[base(data[i][0], j) for j in range(M)] for i in range(N)]
    phi = np.matrix(phi)
    t = np.matrix([x[1] for x in data]).T
    S_N = np.linalg.inv(ALPHA * np.eye(9) + BETA * phi.T.dot(phi))
    m_N = BETA * S_N.dot(phi.T).dot(t)
    
    def f(v):
        s = 0
        for i in range(M):
            s += base(v,i) * m_N[i,0]
        return s
    
    def sigma_N(v):
        phi_x = np.array([])
        for i in range(M):
            phi_x = np.append(phi_x, [base(v,i)])
        return math.sqrt(1/BETA + phi_x.T.dot(S_N).dot(phi_x))

    x = np.arange(0, 1, 0.01)
    y = f(x)
    yt = np.array(list(map(sinpi, x)))
    plt.plot(x, y)
    plt.plot(x, yt)
    sigma = np.array(list(map(sigma_N, x)))
    y_h = y + sigma
    y_l = y - sigma
    plt.fill_between(x, y_h, y_l, facecolor='y', alpha=0.5)
    plt.scatter([x[0] for x in data], [x[1] for x in data], c="red")
    plt.show()

if __name__=='__main__':
    main()