昨天跟大家介紹了貝氏回歸與預測分佈，那我們今天就來看一下程式碼會長什麼樣子吧！
這邊我所使用的資料集是這個

import numpy as np
import matplotlib.pyplot as plt
import scipy.io as io
from numpy.random import multivariate_normal as normal
from math import exp

x = io.loadmat('2_data.mat')['x']
t = io.loadmat('2_data.mat')['t']

def phi(x):
    M = 7
    def sigmoid(x,j):
        s = 0.1
        muj = (2*j)/float(M)
        return 1./ (1. + exp(-1. * ((x - muj)/s)))
    X = []
    for i in range(M):
        X += [sigmoid(x,i)]
    return np.asarray(X)

for n in [10,15,30,80]: #分別用不同的資料量測試
    plt.figure()
    plt.plot(x[:n],t[:n],'o') #劃上資料點

    PHI = phi(x[0])
    target = t[:n]
    for i in range(1,n):
        PHI = np.vstack((PHI, phi(x[i])))

下面這段就是在計算

    s0_inv = (10**-6) * np.identity(7) #S0
    sn_inv = s0_inv + PHI.T.dot(PHI) # 計算Sn，且假設beta = 1
    sn = np.linalg.inv(sn_inv)
    mn = sn.dot(PHI.T.dot(target))

得到 w 的平均與變異數，便可利用numpy提供的normal，隨機產生一些 w 來作圖。

    mn_ = np.reshape(mn, 7)
    line = np.linspace(0., 2., 50)
    ws = normal(mn_, sn, 5)
    for w in ws:
        value = []
        for point in line:
            value += [phi(point).T.dot(w)]
        plt.plot(line, value, linestyle ='--', zorder = 1)
    plt.savefig('%d.png'%n)

接著計算預測值分佈的平均值與變異數，以及畫圖。

    plt.figure()
    plt.plot(x[:n],t[:n],'o')
    mx = []
    vx = []
    for point in line:
        mx += [mn.T.dot(phi(point))] #預測分佈的平均
        vx += [1. + (phi(point).T.dot(sn)).dot(phi(point))] #預測分佈的變異數
    mx = np.reshape(np.asarray(mx), len(mx))
    vx = np.reshape(np.asarray(vx), len(vx))
    plt.plot(line, mx, linestyle = '-', zorder = 1, color = 'red')
    plt.fill_between(line, mx-vx, mx+vx, color = 'pink') #範圍內塗色
    plt.savefig('pred_%d.png'%n)

簡單的實做貝氏回歸後，明天我們就開始跟大家介紹線性模型的另一種問題 — 分類！