贝叶斯曲线拟合实现

100 阅读 0 评论 66 点赞

我是靠谱客的博主寂寞白猫，这篇文章主要介绍贝叶斯曲线拟合实现，现在分享给大家，希望可以做个参考。

贝叶斯预测分布公式： $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)^tssum^{n}_{n=1}phi(x_n)t_n$

$s(x)=beta^{-1}+phi(x)^tsphi(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()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
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()