PRML-1.2.6 贝叶斯曲线拟合

一些记号和回顾

参数	含义
\(N\)	样本量
\(x=(x_1,...,x_N)^T\)	\(样本数据集\)
\(t=(t_1,...,t_N)^T\)	\(样本的目标数据集\)
\(p(x|\mu,\sigma^2)=\prod\limits_{n=1}^{N} \mathcal{N}(x_n|\mu,\sigma^2)\)	\(数据集x是独立同分布，给定\mu和\sigma^2的情况下的数据集的概率\)
\(w\)	\(模型参数\)
\(\mu\)	\(似然函数期望\)
\(\sigma^2\)	\(似然函数方差\)
\(\beta\)	\(似然函数精度,\beta^{-1}=\sigma^2\)
\(y(x,w)=w_0+w_1x+w_2x^2+...+w_Mx^M=\sum\limits_{j=0}^{M}w_jx^j\)	\(多项式拟合函数\)
\(\alpha\)	\(先验分布的精度\)

---

虽然我们已经谈到了先验分布\(p(w | \alpha)\)，但是我们⽬前仍然在进⾏\(w\)的点估计(最大似然估计)，这并不是贝叶斯观点,在⼀个纯粹的贝叶斯⽅法中，我们应该⾃始⾄终地应⽤概率的加和规则和乘积规则。我们稍后会看到，这需要对所有w值进⾏积分。对于模式识别来说，这种积分是贝叶斯⽅法的核⼼

1.纯贝叶斯派的推导和计算

\(简单地说，贝叶斯⽅法就是⾃始⾄终地使⽤概率的加和规则和乘积规则。因此预测概率可以 写成下⾯的形式(假设参数\alpha和\beta是固定的)\)
\(p(t|\mathbb{x},x,t)=\int p(t|x,w)p(w|x,t)dw,-----1.68,\mathbb{x}是新的观测值,要进行预测,x是样本\)
\(p(t|x,w)可以看上一章的1.60公式,p(t|x,w,\beta)=\mathcal{N}(t|y(x,w),\beta^{-1})给出\)
\(p(w|x,t)是参数的后验分布，可以通过公式1.66 \color{red}{p(w|x,t,\alpha,\beta)\propto p(t|x,w,\beta)p(w|\alpha)}归一化得到,\color{red}{后面会学习到这个后验分布是一个高斯分布,可以解系的求出,类似的1.68中积分也可以解系的求解}\)
\(因此，预测分布由⾼斯的形式给出\)
\(p(t|\mathbb{x},x,t)=\mathcal{N}{t|m(x),s^2(x)}\)
\(其中均值和方差分别为\)
\(m(x)=\beta\phi(x)^TS\sum_{n=1}^{N}\phi(x_n)t_n\)
\(s^2(x)=\beta^{-1}+\phi(x)^TS\phi(x)\)
\(\color{red}{这部分没看懂.......不知道怎么推导}\)

2.贝叶斯回归代码

点击查看代码

import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures


# %matplotlib inline
# 定义数据集：
def uniform(size):
    x = np.linspace(0, 1, size)
    return x.reshape(-1, 1)


def create_data(size):
    x = uniform(size)
    np.random.seed(42)
    y = sin_fun(x) + np.random.normal(scale=0.15, size=x.shape)
    return x, y


def sin_fun(x):
    return np.sin(2 * np.pi * x)


x_train, y_train = create_data(10)
x_test = uniform(100)
y_test = sin_fun(x_test)

poly = PolynomialFeatures(9)
X_train = poly.fit_transform(x_train)
X_test = poly.fit_transform(x_test)


#### 定义贝叶斯估计方法：


class BayesianRegressor():

    def __init__(self, alpha=1., beta=1.):
        self.alpha = alpha
        self.beta = beta
        self.mean_prev = None
        self.S = None

    def fit(self, X, t):
        S_inv = self.alpha * np.eye(np.size(X, 1)) + self.beta * np.matmul(X.T, X)
        mean_prev = np.linalg.solve(
            S_inv,
            self.beta * np.matmul(X.T, t)
        )
        self.mean_prev = mean_prev
        self.S = np.linalg.inv(S_inv)

    def predict(self, X):
        y = np.matmul(X, self.mean_prev)
        y_var = 1 / self.beta + np.sum(np.matmul(X, self.S) * X, axis=1)
        y_std = np.sqrt(y_var)
        return y, y_std


# 最后使用我们的模型拟合观察数据，然后对新数据进行预测并绘制曲线：

model = BayesianRegressor(alpha=2e-3, beta=2)
y_train = y_train.reshape(10)
model.fit(X_train, y_train)
y, y_std = model.predict(X_test)
fig = plt.figure(figsize=(12, 8))
plt.scatter(x_train, y_train, facecolor="none", edgecolor="b", s=50, label="training data")
plt.plot(x_test, y_test, c="g", label="$\sin(2\pi x)$")
plt.plot(x_test, y, c="r", label="Mean")
plt.fill_between(x_test[:, 0], y - y_std, y + y_std, color="pink", label="std", alpha=0.5)
plt.title("M=9")
plt.legend(loc=2)
plt.show()