贝叶斯曲线拟合实现

最新推荐文章于 2023-05-16 22:32:47 发布

萌新待开发

最新推荐文章于 2023-05-16 22:32:47 发布

阅读量1.1k

点赞数 2

分类专栏： ᕦ 机器学习 ᕤ 文章标签：概率论机器学习 r语言

本文链接：https://blog.csdn.net/qq_44785318/article/details/120689990

版权

ᕦ 机器学习 ᕤ 专栏收录该内容

68 篇文章 36 订阅

订阅专栏

贝叶斯预测分布公式： $p(y|x,X,Y)=\int p(y|x,w)p(w|X,Y)dw$

这里我们假设后验分布是高斯分布，那么 $p(y|x,X,Y)=\textit{N}(y|m(x),s^2(x))$

其中方差和平均值公式为：

$m(x)=\beta \phi(x)^TS\sum^{N}_{n=1}\phi(x_n)t_n$

$s(x)=\beta^{-1}+\phi(x)^TS\phi(x)$

$S^{-1}=\alpha I+\beta \sum^N_{n=1}\phi(x_n)\phi(x)^T$

import numpy as np
import matplotlib.pyplot as plt


# Number of training points
N = 200
# Precision of targets
beta = 8
# Precision of prior weights distribution
alpha = 0.00

def trueY(X):
    return np.sin(X) * 0.4 +np.cos(X)**2 + (np.cos(X)+1)**(5/6)

def generate(num):
    X = np.random.rand(num)*10
    X = np.sort(np.asarray(X))
    Y = trueY(X)+ np.random.normal(0, beta **-1, size=X.shape)
    return X.tolist(),Y.tolist()

training,y = generate(N)


'''========================================'''
class Bayesian_curve(object):
    def __init__(self,poly = 9):
        self.poly = poly

    def phi(self,data):
        phi_all = []
        for j in data:
            phi_all.append([j ** i for i in range(self.poly + 1)])

        phi_all = np.array(phi_all).T
        I = np.identity(self.poly + 1)
        res = np.sum([i.reshape(-1, 1).dot(i.reshape(-1, 1).T) for i in phi_all.T], axis=0)
        S = np.linalg.inv(beta * res + alpha * I)
        return phi_all, S

    # Mean of Gaussian distribution of t corresponding to x
    def m(self,data,y):
        Phi_all, S = self.phi(data)
        result = np.sum([y[i] * Phi_all.T[i] for i in range(len(y))], axis=0)

        return beta * (Phi_all.T).dot(S.dot(result)).reshape(-1)

    # Variance of Gaussian distribution of t corresponding to x
    def var(self,data):
        Phi_all, S = self.phi(data)
        return np.diagonal((beta ** -1 + (Phi_all.T).dot(S).dot(Phi_all)))
for i in range(5,10):
    bayesian = Bayesian_curve(poly=i)
    plt.fill_between(training, bayesian.m(training,y) - bayesian.var(training),bayesian.m(training,y) + bayesian.var(training), color="g", alpha=0.2)
    plt.plot(training, bayesian.m(training,y), color='r', label='predict')
    plt.scatter(training, y, color='g', alpha=0.3, label='points')
    plt.plot(training, trueY(training), color='b', label="true")
    plt.legend()
    plt.show()