PRML-第1章 绪论 总结

一些数学符号

N
样本量


xn
每一个数据点,或者叫样本点,工程中/训练集中的x_train


tn
训练集中的y_train


CDF
累计概率分布


PDF
概率密度函数


D
观测


X=(X1,...,XN)T

N,x_train


T=(T1,...,TN)T

N,y_train


x=(x1,...,xD)T

1,D,x_test


t=(t1,...,tD)T

1,y_test

w
自然参数

1.1 多项式拟合

y(x,w)=w0+w1x+w2x2+...+wMxM=j=0Mwjxj(1.1)


E(w)=12n=1N{y(xn,w)tn}2(1.2)
误差函数 Loss Function,此处是平方和误差


ERMS=2E(w)/N (1.3)
均方根误差,为了比较不同大小的数据集和保证t有相同单位。引入均方根误差


E~(w)=12n=1N{y(xn,w)tn}2+λ2w2$\tag{1.4}
正则化

w2≡=w02+w12+...+wM2λ的大小控制的正则化影响的大小


w=(XTX)1XTY,w=(XTX)1XTt
最小二乘法解析解


w=(XTX+λI)1XTt
最小二乘法正则化后的解析解,(XTX+λI)w=XTt,

多项式拟合的python实现

点击查看代码

# preparation
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

from prml.preprocess import PolynomialFeature
from prml.linear import (
    LinearRegression,
    RidgeRegression,
    BayesianRegression
)

np.random.seed(1234)


def create_toy_data(func, sample_size, std):
    '''
    产生包含随机高斯噪声的的样本
    :param func:
    :param sample_size:
    :param std:
    :return: x:
    '''
    x = np.linspace(0, 1, sample_size)
    t = func(x) + np.random.normal(scale=std, size=x.shape)  # 产生目标函数的采样,添加高斯分布的噪声,loc 均值,std标准差
    # 这个np.random.normal 产生的随机数每次都一样
    return x, t


def func(x):
    '''
    目标拟合函数
    :param x:
    :return:
    '''
    return np.sin(2 * np.pi * x)


x_train, y_train = create_toy_data(func, 10, 0.25)
x_test = np.linspace(0, 1, 100)
y_test = func(x_test)

# 原目标函数展示
# 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.legend()
# plt.show()

print("x_train =", x_train)

# 几个多项式拟合
for i, degree in enumerate([0, 1, 3, 9]):
    print(i, degree)
    plt.subplot(2, 2, i + 1)
    feature = PolynomialFeature(degree)
    X_train = feature.transform(x_train)  # 输入样本x 做多项式特征处理
    print("X_train.shape=", X_train.shape)
    X_test = feature.transform(x_test)
    print("X_test.shape=", X_test.shape)

    model = LinearRegression()
    model.fit(X_train, y_train)
    y, y_std = model.predict(X_test, True)
    print("y=", y)
    print("y_std", y_std)  # y_std 是X_train 算出来的,和X_test无关

    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="fitting")
    plt.ylim(-1.5, 1.5)
    plt.annotate("M={}".format(degree), xy=(0.8, 1))
plt.legend(bbox_to_anchor=(1.05, 0.64), loc=0, borderaxespad=0.)
plt.show()


点击查看代码

import itertools
import functools
import numpy as np


class PolynomialFeature(object):
    """
    多项式组合

    将输入的数组变成多项式组合

    Example
    =======
    x =
    [[a, b],
    [c, d]]

    y = PolynomialFeatures(degree=2).transform(x)
    y =
    [[1, a, b, a^2, a * b, b^2],
    [1, c, d, c^2, c * d, d^2]]
    """

    def __init__(self, degree=2):
        """
        构建多项式特征

        参数
        ----------
        degree : int
            多项式的阶数
        """
        assert isinstance(degree, int)  # 判断 degree 是不是int类型
        self.degree = degree

    def transform(self, x):
        """
        将输入的数组变成多项式组合

        参数
        ----------
        x : (sample_size, n) ndarray
            input array

        Returns
        -------
        output : (sample_size, 1 + nC1 + ... + nCd) ndarray
            polynomial features
        """
        if x.ndim == 1:  # ndim 矩阵x的维度,是几维矩阵
            x = x[:, None]
        x_t = x.transpose()  # 转置
        features = [np.ones(len(x))]  # 以x的形状做全1 矩阵
        for degree in range(1, self.degree + 1):
            for items in itertools.combinations_with_replacement(x_t,
                                                                 degree):  # 排列组合 这方法允许 元素中的数据重复出现 https://blog.csdn.net/weixin_43866211/article/details/101758931
                features.append(functools.reduce(lambda x, y: x * y, items))
        return np.asarray(features).transpose()


点击查看代码

import numpy as np
from prml.linear.regression import Regression


class LinearRegression(Regression):
    """
    线性回归模型
    y = X @ w - 模型公式
    t ~ N(t|X @ w, var) - 噪声服从正态分布
    """

    def fit(self, X: np.ndarray, t: np.ndarray):
        """
        执行最小二乘拟合

        参数
        ----------
        X : (N, D) np.ndarray
            训练集的观测数据
        t : (N,) np.ndarray
            训练集的目标数据
        """
        # 直接得到w的解析解
        self.w = np.linalg.pinv(X) @ t  # @ 运算是矩阵相乘,若是向量就是內积,见下面的main函数
        # np.linalg.pinv 是求伪逆,公式为 pinv(X) = (X^TX)^{−1}X^T
        self.var = np.mean(np.square(X @ self.w - t))  # 均方误差 (X^Tw -t)^2

    def predict(self, X: np.ndarray, return_std: bool = False):
        """
        Given input的情况下做预测

        Parameters
        ----------
        X : (N, D) np.ndarray
            预测的输入 是一组输入,N是输入的数据数量,D是多项式的阶数
        return_std : bool, optional
            是否需要返回均方误差

        Returns
        -------
        y : (N,) np.ndarray
            预测值
        y_std : (N,) number
            均方误差
        """
        y = X @ self.w
        if return_std:
            # y_std = np.sqrt(self.var) + np.zeros_like(y)
            y_std = np.sqrt(self.var)  # 这里的var是X_train 算出来的均方误差
            return y, y_std
        return y


if __name__ == '__main__':
    X = np.array([1, 2, 3, 4])
    w = np.array([6, 7, 8, 9])
    print(X @ w)

    X = np.array([[1, 2], [1, 2]])
    w = np.array([[6, 7], [6, 7]])
    print(X @ w)


点击查看代码

class Regression(object):
    """
    Base class for regressors
    """
    pass


均方根误差 python

点击查看代码


def rmse(a, b):
    return np.sqrt(np.mean(np.square(a - b)))


training_errors = []
test_errors = []

for i in range(10):
    feature = PolynomialFeature(i)
    X_train = feature.transform(x_train)
    X_test = feature.transform(x_test)

    model = LinearRegression()
    model.fit(X_train, y_train)
    y = model.predict(X_test)
    training_errors.append(rmse(model.predict(X_train), y_train))
    test_errors.append(rmse(model.predict(X_test), y_test + np.random.normal(scale=0.25, size=len(y_test))))

plt.plot(training_errors, 'o-', mfc="none", mec="b", ms=10, c="b", label="Training")
plt.plot(test_errors, 'o-', mfc="none", mec="r", ms=10, c="r", label="Test")
plt.legend()
plt.xlabel("degree")
plt.ylabel("RMSE")
plt.show()

正则化python实现

点击查看代码

feature = PolynomialFeature(9)
X_train = feature.transform(x_train)
X_test = feature.transform(x_test)

model = RidgeRegression(alpha=1e-3)
model.fit(X_train, y_train)
y = model.predict(X_test)

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="fitting")
plt.ylim(-1.5, 1.5)
plt.legend()
plt.annotate("M=9", xy=(-0.15, 1))
plt.show()

1.2概率论

nij
N次取样,每次得到一组的x,y。X=xi,Y=yi的取到次数是nij


(1.5)p(X=xi,Y=yi)=nijN (1.6)p(X=xi)=ciN
基本公式


p(X=xi)
边缘概率


(1.8)p(Y=yj|X=xi)=nijci
条件概率


(1.7)p(X=xi)=j=1Lp(X=xi,Y=yi)(1.10)p(X)=Yp(X,Y)(1.31)p(X)=Yp(X,Y)dy


|(1.9)p(X=xi,Y=yj)=nijN=nijciciN=p(Y=yj|X=xi)p(X=xi)(1.11)p(X,Y)=P(YX)p(X)
乘积法则


(1.12)P(Y|X)=P(XY)P(Y)P(X)(1.12)P(Y|X)=P(XY)P(Y)YP(X|Y)P(Y)(1.12)P(Y|X)=P(X|Y)P(Y)YP(X|Y)P(Y)dy
贝叶斯定理


P(X)=YP(X|Y)P(Y)P(X)=YP(X|Y)P(Y)dy
归一化常数


P(Y)
先验概率


P(Y|X)
后验概率


(1.27)py(y)=px(x)||dxdy=px(g(y))g(y)
Jacobian因子变化(概率形式变换),有点像复合函数,只不过这个函数是个概率密度函数,推导见其他文章


(1.33)E[f]=xp(x)f(x)
离散函数的期望


(1.34)E[f]=p(x)f(x)dx
连续函数的期望

(1.35)E[f]1Nn=1Nf(xn)
大数定理

(1.36)Ex[f(x,y)]
多变量下对x变量求期望


(1.37)E[f|y]xp(x|y)f(x)
条件期望


(1.38)var[f]=E[(f(x)E[f(x)])2](1.39)var[f]=E[f(x)2]E[f(x)]2(1.40)var[x]=E[x2]E[x]2
方差,第三个概率论书上都有证明


(1.41)cov[x,y]=Ex,y[xE[x]yE[y]] =Ex,y[xy]E[x]E[y](1.42)cov[x,y]=Ex,y[xE[x]yTE[yT]] =Ex,y[xyT]E[x]E[yT]
协方差,第二个是向量表示

1.23 贝叶斯概率

p(D|w)
似然函数

p(w|D)
后验

p(w)
先验

(1.45)p(D)=p(D|w)P(w)dwp(D)=wp(D|w)P(w)
归一化因子

(1.43)p(w|D)=p(D|w)p(w)p(D)
贝叶斯定理

(1.44)posteriorlikelihood×prior
三者关系

频率派:w是固定的参数。最小误差函数就是最大似然估计


贝叶斯派:w不是固定的,需要用概率分布表达这种不确定性。

1.24 高斯分布

(1.46)N(x|μ,σ2)=1(2πσ2)12exp{12σ2(xμ)2}
一元变量的高斯分布


(1.52)N(x|μ,Σ)=1(2π)D/21Σ1/2exp12(xμ)TΣ1(xμ)
D维高斯分布


Σ
D×D维协方差矩阵


(1.49)E[x]=N(x|μ,σ2)xdx=μ
高斯分布的均值-一阶矩


(1.50)E[x2]=N(x|μ,σ2)x2dx=μ2+σ2
二阶矩


(1.51)var[x]=E[x2]E[x]2=σ2
方差


样本均值:(1.55)μML=1Nn=1Nxn
样本方差:(1.56)σML2=1Nn=1N(xnμML)2
高斯分布的最大似然函数最优解

--
(1.57)E[μML]=μ
高斯分布的均值\mu是无偏估计


(1.58)E[σML2]=(N1N)σ2
高斯分布的方差是有偏估计


(1.59)σ~2=NN1σML2=1N1n=1N(xnμML)2
这个方差参数估计是无偏的


Var(X1+X2+...+Xn)=k=1nVar(Xk)
X1,X2,...,Xn

1.25 重新考察曲线拟合-从后验(MAP)角度看似然问题

什么是后验?就是从验证数据集 x,t角度看最大似然估计

X=(X1,...,XN)T

N,x_train


T=(T1,...,TN)T

N,y_train


x=(x1,...,xD)T

1,D,x_test


t=(t1,...,tD)T

1,y_test


β
精度,β1=σ2


(1.60)p(tx,w,β)=N(ty(x,w),β1)

y,tp(tx,w,β),x,t


(1.61)p(TX,w,β)=n=1NN(tny(xn,w),β1)

似然函数


(1.62)lnp(TX,w,β)=β2n=1N{y(xn,w)tn}2+N2lnβN2ln(2π)

对式1.60进行最大似然,发现最大化该似然函数等价于最小化误差函数


(1.63)1βML=1Nn=1Ny(xn,wML)tn2

最大似然解


(1.64)p(t|x,wML,βML)=N(t|y(x,wML),βML1)

代入最大似然解后的预测分布


(1.65)p(w|α)=N(w|0,α1I)=(α2π)M+12exp{α2wTw}
w的先验,假设是正态的(1.69中可以证明先验是正态的)


α
w,


(1.66)p(w|X,T,α,β)p(T|X,w,β)p(w|α)
证明见本博客


(1.67)β2n=1Ny(xn,w)tn2+α2wTw
对1.66取log得到MAP公式,就是正则化函数

1.26 纯贝叶斯方法的曲线拟合


(1.68)p(t|x,X,T)=p(t|x,w)p(w|X,T)dw
本章尝试对w积分,而不是像上一章对w做估计(或者说是仅利用加法法则和乘法法则求后验,而不是用贝叶斯定理求后验)


1.68求解结果
(1.69)p(t|x,X,T)=N(t|m(x),s2(x))
(1.70)m(x)=βϕ(x)TSn=1Nϕ(xn)tn
(1.71)s2(x)=β1+ϕ(x)TSϕ(x)
(1.72)S1=αI+βn=1Nϕ(x)ϕ(x)T
其中I是单位矩阵,定义向量ϕ(x)ϕi(x)=xi,i=0,...,M

总结一下几个正态分布
1.预测分布-MAP
(1.64)p(t|x,wML,βML)=N(t|y(x,wML),βML1)
2.先验分布-MAP
(1.65)p(w|α)=N(w|0,α1I)=(α2π)M+12exp{α2wTw}
3.先验分布
(1.68中的推导可得)p(w|X,T)N(w|μ,S)
S1=αI+βϕTϕ,μ=βSϕTtn
4.预测分布
(1.69)p(t|x,X,T)=N(t|m(x),s2(x))

贝叶斯回归的python实现

注意,在这段python代码中,实际没有用到上面的1.68-1.72的公式,实际是共轭先验的思想,对w假设是高斯分布,则其对应的后验也是高斯,推导的结果见公式3.48-3.51,最后样本拟合后得到的是p(w|X,T)这样一个正态分布,注意是一个分布,然后针对这个分布做随机采样,就能得到y_test

(3.48)p(w)=N(w|m0,S0)
(3.49)p(w|t)=N(w|mN,SN)
其中
(3.50)mN=SN(S01m0+βΦTt)(3.51)SN1=S01+βΦTΦ

版本一

点击查看代码
import numpy as np
from prml.linear.regression import Regression


class BayesianRegression(Regression):
    """
    贝叶斯回归模型

    w ~ N(w|0, alpha^(-1)I)
    y = X @ w
    t ~ N(t|X @ w, beta^(-1))
    """

    def __init__(self, alpha: float = 1., beta: float = 1.):
        self.alpha = alpha
        self.beta = beta
        self.w_mean = None
        self.w_precision = None

    def _is_prior_defined(self) -> bool:
        return self.w_mean is not None and self.w_precision is not None

    def _get_prior(self, ndim: int) -> tuple:  # -> tuple 表示返回一个tuple
        '''
        :param ndim:
        :return:  w_mean  先验的均值
                    w_precision 先验的精度
        '''
        if self._is_prior_defined():
            return self.w_mean, self.w_precision
        else:
            return np.zeros(ndim), self.alpha * np.eye(ndim)

    def fit(self, X: np.ndarray, t: np.ndarray):
        """
        给定训练集,更新贝叶斯参数

        Parameters
        ----------
        X : (N, n_features) np.ndarray
            训练集输入
        t : (N,) np.ndarray
            训练集输出
        """

        # 先验的均值,先验的精度
        mean_prev, precision_prev = self._get_prior(
            np.size(X, 1))  # mean_prev =m_0 precision_prev =S_0= \alpha I 见PRML 3.48

        w_precision = precision_prev + self.beta * X.T @ X  # 这里是计算 w_precision=S_N^{-1}=S_0^{-1}+\beta X^T X,PRML 3.54

        # 计算m(x) 也就是后验的均值
        w_mean = np.linalg.solve(
            w_precision,
            precision_prev @ mean_prev + self.beta * X.T @ t
            # w_mean=m_N ,解方程 S_N^{-1} w_mean =S_0^{-1} * m_0 +\beta X^T t ,PRML 3.53
        )
        self.w_mean = w_mean
        self.w_precision = w_precision
        self.w_cov = np.linalg.inv(self.w_precision)  # w_cov = S_N , 协方差矩阵

    def predict(self, X: np.ndarray, return_std: bool = False, sample_size: int = None):
        """
        return mean (and standard deviation) of predictive distribution

        Parameters
        ----------
        X : (N, n_features) np.ndarray
            independent variable
        return_std : bool, optional
            flag to return standard deviation (the default is False)
        sample_size : int, optional
            number of samples to draw from the predictive distribution
            (the default is None, no sampling from the distribution)

        Returns
        -------
        y : (N,) np.ndarray
            mean of the predictive distribution
        y_std : (N,) np.ndarray
            standard deviation of the predictive distribution
        y_sample : (N, sample_size) np.ndarray
            samples from the predictive distribution
        """

        # 得到拟合的正态分布后,对这个正态分布采样
        if sample_size is not None:
            w_sample = np.random.multivariate_normal(
                self.w_mean, self.w_cov, size=sample_size
            )
            y_sample = X @ w_sample.T
            return y_sample
        y = X @ self.w_mean
        if return_std:
            y_var = 1 / self.beta + np.sum(X @ self.w_cov * X, axis=1)
            y_std = np.sqrt(y_var)
            return y, y_std
        return y

点击查看代码

# 贝叶斯曲线拟合
model = BayesianRegression(alpha=2e-3, beta=2)
model.fit(X_train, y_train)

y, y_err = model.predict(X_test, return_std=True)
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, y - y_err, y + y_err, color="pink", label="std.", alpha=0.5)
plt.xlim(-0.1, 1.1)
plt.ylim(-1.5, 1.5)
plt.annotate("M=9", xy=(0.8, 1))
plt.legend(bbox_to_anchor=(1.05, 1.), loc=2, borderaxespad=0.)
plt.show()


版本二

点击查看代码

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()


1.5.4 推断和决策

三种方法
1.生成模型:通过对每个类别Ck,独立的确定类别的条件密度p(x|Ck)来解决推断问题,还分别推断出类别的先验概率p(Ck),然后使用贝叶斯定理:

(1.82)p(Ck|x)=p(x|Ck)p(Ck)p(x)

来计算类别的后验概率p(Ck|x)
2.判别模型,解决确定类别的后验密度p(Ck|x)的推断问题,然后,使用决策论来对新的输入x进行分类。
3.判别函数:找到能直接把输入x映射到类别标签f(x)

1.6 信息论

离散熵:如果箱子是X的离散状态xi,X的熵就是:

(1.98)H[p]=ip(xi)lnp(xi)

当所有的p(xi)都相等,且值为p(x_i)=\frac{1}{M}时,熵取得最大值
微分熵

(1.104)H[x]=p(x)lnp(x)dx

最大化微分熵的分布是高斯分布
高斯分布的微分熵

(1.110)H[x]=12{1+ln(2πσ2)}

条件熵

(1.111)H[y|x]=p(y,x)lnp(y|x)dydx

条件熵的乘积规则

(1.112)H[x,y]=H[y|x]+H[x]

KL散度(相对熵),KL散度不是对称的,KL(p|q)KL(q|p)

(1.113)KL(pq)=p(x)lnq(x)dx(p(x)lnp(x)dx)=p(x)ln{q(x)p(x)}dx

互信息

(1.120)I[x,y]KL(p(x,y)p(x)p(y)) =p(x,y)ln(p(x)p(y)p(x,y))dxdy

利用概率加和和乘积规则

(1.121)I[x,y]=H[x]H[x|y]=H[y]H[y|x]

posted @   筷点雪糕侠  阅读(147)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)
· spring官宣接入deepseek,真的太香了~
点击右上角即可分享
微信分享提示