[1]崔佳旭,杨博.贝叶斯优化方法和应用综述.软件学报,2018,29(10):3068-3090. http://www.jos.org.cn/1000-9825/5607.htm
[2]Bobak Shahriari, Kevin Swersky,Ziyu Wang, Ryan PAdams, and Nando de Freitas. 2016. Taking the human out of the loop: A review of bayesian optimization. Proc. IEEE 104,1(2016),148–175.
简介
首先，贝叶斯优化能干什么？给我的感觉是无所不能，当然其效果有些可能不尽如人意。贝叶斯优化，可以做回归的东西（虽然总感觉这个东西只是一个附属品），然而主要是去解决一个“优化问题”。
贝叶斯优化解决的是下面类型的问题：
注: 使用"argmin"并无实质上的不同，事实上[1]中采用的便是"argmin"。
往往，(f)我们并不知道，所以，这类问题很难采用经典的梯度上升("argmin"则梯度下降)来解决。贝叶斯优化采用概率代理模型来应对。(mathrm{x})是决策，往往称(mathcal{X})为决策空间。药物配方是一种决策，神经网络卷积核大小等也可以看成一种决策。而且，这种决策与最后的输出的关系（即(f)）往往很难知晓。这也正是贝叶斯优化的强大之处。
贝叶斯优化框架
上面俩幅图分别来自[2]和[1]，因为一些符号的差异，往下除特别指明，采用的均为[2]中的符号。
贝叶斯优化，每一次迭代，首先在代理模型的“先验”下，通过最大化采集函数（该函数往往是对评估点的分布以及(f(x))的提升的一种权衡(trade-off)）。新的评估点，作为输入传入系统，获得新的输出，以此来更新(D)和概率代理模型。
其中(D_t={(mathrm{x_1}, y_1), ldots, (mathrm{x_t},y_t)})
上面这幅图，便是贝叶斯优化的一个简单演示。黑色虚线表示目标函数(f),而黑色实线表示我们拟合的曲线（图中是通过对概率代理模型求均值获得的）。蓝紫色区域是(mu(cdot)pm sigma(cdot))。下方的绿色曲线则是每一次迭代的(alpha(cdot))，可以看出，每一次迭代选出的评估点都是(alpha(cdot))最大值所对应的(mathrm{x})。
下面，我们分别来讨论概率代理模型，以及采集函数。
概率代理模型
概率代理模型，顾名思义，就是用来代理(f(cdot))的一个概率模型。
参数模型
参数模型，即(f(cdot))可由参数(mathrm{w})来决定。如果我们给定(mathrm{w})的先验分布(p(mathrm{w}))。那么，通过贝叶斯公式，我们可以获得(mathrm{w})的后验分布：
现在问题来呢，我们还不知道(p(D|mathrm{w}))和(p(D))啊。(p(D|w))是一个似然分布，往往通过(prod p((mathrm{x}_i,y_i)|mathrm{w}))来计算，当然，我们得知道(p((mathrm{x}_i, i_i)|mathrm{w}))。至于(p(D))，比较难计算，但是，(p(D))在这里只是扮演了系数的作用，所以用核方法就能解决。事实上，我们常常选择共轭先验分布作为(mathrm{w})的先验分布。
汤普森采样和Beta-Bernouli模型
这里给出一个例子：实验室有K种药，我们需要通过药物实验来找出哪种药的效果最好。假设，药作用在某个病人身上只有成功治愈和失败俩种可能，且不能通过一种药的效果来判断另外一种药的疗效。这种类型的问题似乎被称为A/B测试，常用于新闻推荐等。
我们用(a in 1 ldots K)来表示药物，(w_a)表示第(a)种药物成功治愈病患的可能性，而(y_i in {0, 1})则表示病人(i)的治疗情况（0失败，1治愈）。函数(f(cdot))就是某种复杂的映射。让参数(mathrm{w} in (0, 1)^{K})，(mathrm{w}_a=w_a)。那么我们选择的概率代理模型是(f_mathrm{w}(a):=w_a)。
我们选择(mathrm{Beta})分布作为(mathrm{w})的先验分布，因为这是其共轭先验分布。
定义：
其中(n_{a,0})表示(n)次评估中，选中(a)药物且治疗失败的数目，(n_{a,1})则反之。(mathbb{I}(cdot))只有((cdot))成立为1否则为0。
那么，(mathrm{w})的后验概率为：
上述推导见附录。
从上述也能发现，超参数(alpha, eta)表示的治愈数和失败数。下图是以(mathrm{Beta}(w|2,2))为先验的一个例子。
汤普森采样-wiki
那么在(D_n)的基础上，如果找(a_{n+1})呢。以下采用的是汤普森采样（或后验采样）：
(widetilde{mathrm{w}}sim p(mathrm{w}|D_n)),即(widetilde{mathrm{w}})从(mathrm{w})的后验分布中采样得到。
该模型的好处是：
- 只有(alpha,eta)俩个超参数
- 是对exploration和exploitation的一种比较好的权衡
- 比较容易和蒙特卡洛采样-HFUT_qianyang结合
- 汤普森采样的随机性使其容易推广到批量处理（不懂）
下面是该模型的算法：
线性模型(Linear models)
上述的模型在应对组合类型的时候会显得捉襟见肘。比方，我们在从(d)个元素中寻求一种搭配，每个元素有({0, 1})俩种状态，那么总共就有(2^K)种组合，如果为每种组合都设立一个(w)，显然不切实际。更何况，先前模型的假设（无法从一种组合推断另外一种组合的有效性）显得站不住脚。因为，不同组合往往有微妙的相关性。
采用线性模型，能比较好的解决这一问题。
我们把每一种策略设为(mathrm{x}_a in mathbb{R}^{d})，并且概率代理模型为(f_{mathrm{w}}(a)=mathrm{x}_a^{mathrm{T}}mathrm{w})，现在(mathrm{w})成了权重向量。这只是代理模型的一部分，因为并没有体现“概率”的部分。
组合(a)的观测值为(y_a)，服从正态分布。很自然地，我们同样选择共轭先验分布作为(mathrm{w}, sigma^2)的先验分布：normal_inverse_gamma-wiki
(mathrm{NIG})分布有4个超参数，而(mathrm{w}, sigma^2)的后验分布（(D_n)的条件下）满足：
其中(mathrm{X})的第(i)行为(mathrm{x}_{alpha_i})。
推导见附录。
关于(alpha_{n+1})的选择，同样可以采用汤普森采样：
其中(widetilde{mathrm{w}} sim p(mathrm{w}|D_n))
线性模型有很多扩展：
(f(x)=Phi(mathrm{x})^{mathrm{T}}mathrm{w})
(f(x)=g(mathrm{x^{T}w}))
其中，(Phi(mathrm{x})=(phi_1(mathrm{x}),ldots,phi_J(mathrm{x}))^{mathrm{T}})，(phi_j(mathrm{x}))常常为：
和
这里，(Lambda)，({mathrm{z}_j }_{j=1}^{J})，({ w_j}_{j=1}^{J})均为超参数，至于这些超参数怎么更新，我不大清楚。
非参数模型
非参数模型不是指没有参数，而是指参数（数量）不定。
我们先来看如何把先前的线性模型转换成非参数模型。
我们假设(sigma^2)是固定的，且(p(mathrm{w}|V_0)=mathcal{N}(0,V_0)),即服从均值为(0),协方差矩阵为(V_0)的多维正态分布。那么，(ysim mathcal{N}(mathrm{Xw},sigma^2mathrm{I}))，我们可以积分掉(mathrm{w})从而获得(y)的一个边际分布:
推导见附录。
就像先前已经提到过的，我们可以引入(Phi = (phi_i,ldots,phi_J)^{T}),
将资料（设计）矩阵(X)映射到(mathbb{R}^{n imes J})，如此一来，相应的边际分布也需发生变化：
注意到(Phi V_0 Phi^{T})，事实上，我们不需要特别指明(Phi)，而只需通过kernel.
(K)将是(Phi V_0 Phi^{T})的一个替代。采用这个策略，比原先在计算上和可解释性上更有优势。
不过，还有另外一个问题，如果去寻找下一个评估点呢。寻找下一个评估点，需要我们做预测，但是上面的边际分布显然是无法进行预测的。不过，我们可以通过条件公式来获得：
(mathrm{X}_{*})是新的位置,而(y_{*})是相应的预测，二者都可以是向量。
分子部分是一个联合的高斯分布。到此，我们实际上完成了一个简易的高斯过程，下面正式介绍高斯过程。
高斯过程
高斯过程(f(mathrm{x}) sim GP(mu_0,k))，其核心便是均值函数(mu_0)(在贝叶斯优化中，我们常常选择其为0)和协方差函数(k(mathrm{x}_i,mathrm{x}_j))，而观测值(y=f+varepsilon)。通过高斯过程得到的序列(f_{1:n})，以及观测值(y_{1:n})都服从联合正态分布：
其中(m_i := mu_0(mathrm{x}_i),K_{i,j}:=k(mathrm{x}_i, mathrm{x}_j)),(sigma^2)是随机变量(varepsilon)的方差。
于是，我们可以像之前所做的那样，求边际分布，和(y_*)的分布。
首先，
我们并没有给出这个证明，因为这个不难验证。接下来，为了预测，我们需要求后验分布。论文此时并没有选择(y_*)，而是(f_*)的后验分布，这点倒是可以理解，比较我们的目标就是最大化(f(mathrm{x}))。不过，论文给出的是标量（也就是只有一个预测值），实际上，很容易扩展到多个，在附录里，我们给出多个的后验分布的推导。
我们先给出一个的后验分布，依旧是正态分布：
常用的一些kernels
Kernel method - wiki
Matern covariance function - wiki
文献[1]给出的形式如下(实际上是(d=1)的情况)：
其中，(r=|x-x'|)，(v)为平滑参数，(l)为尺度参数，(K_v)为第二类变形贝塞尔函数。
同时给出了几种常用的Matern协方差函数。
文献[2]给出的是另外一种表示方式：
其中，(r^2=mathrm{(x-x')^{T}Lambda(mathrm{x-x'})})，(Lambda)是一个对角矩阵，其对角线元素为( heta_i^2,i=1,ldots,d)。
这些参数可以这么理解：
- (v)，被称为平滑参数，因为，从使用Matern协方差函数的高斯过程中采样的目标函数(f(x))是(lfloor v-1 floor)次均方可微的（为什么？）。
- (l)和( heta)的功能相似，反映该维度的重要性（前者是越小越重要，后者越大越重要）。
上面的一些参数，会在下面给出一些更新的方法。
边际似然
log 边际似然函数可以表示为：
从图中可以看到，等式右边被分成了三部分，三者有不同含义：
- 第一项，表示模型和数据的拟合程度的好坏，以马氏距离为指标；
- 第二项，表示模型的复杂度；
- 第三项，是数据点(n)的线性函数，表示数据的不确定性随着数据增大而减少（？）。
一个非常自然的想法是，对上述似然函数进行极大似然估计，从而获得( heta:=( heta_{0:d},mu_0,sigma^2))的估计。
复杂度
每一次高斯过程的复杂度在(O(n^3))级别左右，这是由计算矩阵的逆所带来的。通过Cholesky分解，可以降为(O(n^2))。
所以产生了一些算法，试图克服这个困难。
sparse pseudo-input Gaussian processes (SPGP)
SPGP从n个输入中选择m个伪输入替代，从而达到降秩的目的。同时，这m个伪输入的位置也作为参数（虽然我不知道怎么去更新）。好处自然是，
能够把复杂度降为(0(nm^2+m^3))。缺点是，参数相对比较多，容易产生“过拟合”现象。
Sparse spectrum Gaussian processes(SSGP)
由Bochner定理得，任何稳定得kernel都有一个正定得傅里叶谱表示：
之后，通过蒙特卡洛方法，采样m个样本频率，来近似估计上诉的积分。从而获得近似的协方差函数，当数据集较小时，SSGP同样易产生“过拟合”现象。
随机森林
随机森林可以作为高斯过程的一种替代。缺点是，数据缺少的地方，估计的并不准确（边际更是常数）。另外，由于随机森林不连续，也就不可微，所以无法采用梯度下降（或上升）的方法来更新参数。另外不解的是，随机森林的参数，即便我们给了一个先验分布，可其后验分布如何求呢？
采集函数
首先我们有一个效用函数(U:mathbb{R}^d imes mathbb{R} imes Theta ightarrow mathbb{R})，顾名思义，效用函数，是反映评估(mathrm{x})和对应的函数值(v=f(mathrm{x}))（在( heta)条件下）的一个指标。论文[1]并没有引入这个效用函数，论文[2]引入这个概念应该是为了更好的说明。
一种采集函数的选择，便是期望效用：
其实蛮奇怪的，因为对(v)求期望也就罢了，这个采集函数对( heta)也求了期望，我的理解是，这样子更加“模糊”了，如果选择极大似然等方式产出的“精准”的( heta)，或许不能够很好的让评估点足够分散，而陷入局部最优。而且，这样子做，似乎就没有必要去估计参数( heta)，虽然代价是求期望。
从下面的一些算法中我们可以看到，往往没有(mathbb{E}_{ heta}(cdot))这一步骤。
最后再次声明，采集函数的设计，往往都是对exploration和exploitation的一种权衡。即，我们希望新的评估点(mathrm{x})既要和原来的那些数据点远一些（对未知区域的探索），又能够让(f(mathrm{x}))能够提升（对当前区域的开发探索）。
基于提升的策略
PI (probability of improvement)
PI的采集函数的设计思想很简单，就是我们要寻找一个评估点(mathrm{x})，这个(mathrm{x})使得(v=f(mathrm{x}))较已知的最大的（如果一开始是argmin就是最小的）(f(mathrm{x}))，令其为( au)。往往，( au=min_{mathrm{xin X}} f(mathrm{x}))。
采集函数为：
注意，这里的(Phi)是标准正态分布的概率函数。
这个采集函数里的效用函数是：
其中(mathbb{I})为指示函数。
当( au)就是(f(mathrm{x}))的最小值时，PI的效果非常好。
PI一个“弊端”是，只在乎提升的概率，而在乎提升的幅度，而，EI就涵盖了这俩方面。
EI(expected improvement)
通常，其提升函数由下式表示：
而相应的的采集函数是：
其中(phi)是标准正态分布的概率密度函数。式子通过积分变量替换可以推得。
实际上(I)就是效用函数(U)。
UCB(upper confidence bound)
采集函数为：
这个采集函数，可以这么理解，对于任意一个(mathrm{x})，它有一个均值(mu_n (mathrm{x}))，有一个标准差(sigma_n(mathrm{x}))（体现浮动范围和程度），(eta_nsigma_n(mathrm{x}))我们认为比较可靠的界，我们认为，(f(mathrm{x}))有较大可能达到(mu_n(mathrm{x})+eta_nsigma_n(mathrm{x}))的值。所以最大化采集函数，就是最大化我们的这一种希望。
论文[2]中说(eta)的选择往往是Chernoff-Hoeffding界。听起来很玄乎，但是，UCB现在貌似非常火。另外，有一套理论，能够引导和规划超参数(eta_n)，使得能够达到最优。
基于信息的策略
不同之前的策略，基于信息的策略，依赖全局最优解(mathrm{x}^*)的后验分布(p_*(mathrm{x}|D_n))。该分布，隐含在函数(f)的后验分布里（不同的(f)有不同的全局最优解，从而(mathrm{x}^*)也有一个后验分布）。
熵搜索算法旨在寻找能够极大减少不确定度的评估点(mathrm{x}^*)，从另外一个角度来说，(mathrm{x}^*)给与了我们最多的信息。
因此效用函数为：
此时的采集函数为：
其中(H(mathrm{x^*|D_n})=int p(mathrm{x^*|D_n})log p(mathrm{x^*}|D_n)dmathrm{x^*})代表微分熵，(y sim mathcal{N}(mu_n(mathrm{x}), sigma^2_n(mathrm{x})+sigma^2))。
不过，估计(p(mathrm{x*}|D_n))或者(H)都绝非易事，计算量也让人望而却步。
一种新的方法PES(predictive entropy search)很好地改善了这些问题。
利用互信息的对称性（不懂），我们可以将(alpha_{ES}(mathrm{x}))改写为：
使得我们不必苦于(p(mathrm{x^*}|D_n))上。
采集函数的组合
不同采集函数的组合或许能够达到更好的特性。
一种方案是，不同采集函数会给出一系列最优评估点的候选。通过一个准则来判断孰优孰劣。
ESP（entropy search portfolio）采用的准则如下：
即选择那个让不确定性最小的候选点。
出于实际的考量
处理超参数
处理（优化）超参数一般有如下的策略：
- 通过(y)的边际分布，进行极大似然估计，获得(hat{ heta}_n^{ML})
- 给参数( heta)一个先验分布(p( heta))，通过贝叶斯公式获得其后验分布(p( heta | D_n))，可得最大后验估计(hat{ heta}_n^{MAP})
- 同样用到后验分布(p( heta|D_n))，对( heta)进行边际化处理:
这部分不必再估计参数( heta)，代价是需要估计一个期望，这项工作有不同的方案：蒙特卡洛技术，贝叶斯蒙特卡洛技术等。
采集函数的优化
实际上，采集函数的优化并不简单，因为其往往非凸。有基于梯度的方法（如果可微的话），还有网格搜索，自适应协方差矩阵进化策略，多启动局部搜索等。
最近还流行用空间切分的方法替代贝叶斯模型，如SOO（simultaneous optimistic optimization）。在有可用的先验知识时，这类方法比不上贝叶斯优化有效，BamSOO(Bayesian multi-scale SOO)算法结合贝叶斯和树切分空间，有很好的效果。
数值实验
Beta-Bernoulli
我们准备了8枚药丸，治疗概率分别0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8。进行80次试验，每次试验治疗3个无差别病患。下图是一个例子：
再下面的例子是，每次试验不选择药丸，统统进行试验：
高斯过程
选择最简单的一维高斯过程，超参数有( heta_0, heta_1, sigma),选择的kernel为 ( heta_0^2 exp{-frac{1}{2}r^2 })，采集函数为(EI)。优化超参数使用梯度下降（因为别的方法不怎么会），优化采集函数使用的是网格搜索。另外，我们给输出附加了方差为0.0001（基本上没有）的白噪声。
采用[1]中的例子：
首先，是在我们不知道( au)的情况下，我们取初始值0.1，0.2， 0.4，0.5，0.7，0.9.，上面的图是选取11个点后的均值函数，下面的是真实的曲线，上面的点是每次选取的点。
接下来，我们固定( au=-0.2)
最后再来一个加入方差为0.0025白噪声的( au=-0.2)和未知的的：
代码
# Beta-Bernouli
import numpy as np
import scipy.stats
class Pill:
"""模拟药丸
>>> p = Pill(0.5, -1, 1)
Traceback (most recent call last):
...
AssertionError: -1 is not positive
>>> p = Pill(0.5, 1, 0.5)
>>> p = Pill(0.7, 2, 2)
>>> p.patient_num
4
>>> p.cure_num
2
>>> p.fail_num
2
"""
def __init__(self, cure_pro, cure_num, fail_num):
assert 0 <= cure_pro <= 1, "Invalid probability"
Pill.positive(cure_num)
Pill.positive(fail_num)
self.__cure_pro = cure_pro #私有变量 不允许修改
self.__patient_num = cure_num + fail_num
self.__cure_num = cure_num #预先给定的治疗成功的数
self.__fail_num = fail_num #预先给定的治疗失败的数
self.guess_pro = 0 #采样概率
self.__expect_pro = self.__cure_num / self.__patient_num #期望概率
@classmethod
def positive(cls, num):
assert num > 0,
"{0} is not positive".format(num)
@property
def expect_pro(self):
"""获取期望概率"""
return self.__expect_pro
@property
def cure_pro(self):
"""获得cure_pro"""
return self.__cure_pro
@property
def patient_num(self):
"""返回治疗的人数"""
return self.__patient_num
@property
def cure_num(self):
"""获得药丸治愈的人数"""
return self.__cure_num
@property
def fail_num(self):
"""返回失败的人数"""
return self.__fail_num
def curing(self, number=1):
"""使用该药丸治疗病人
返回(成功数，失败数)
"""
if number < 0:
raise ValueError("Invalid number")
self.__patient_num += number
fail = 0
for i in range(number):
if np.random.rand() > self.cure_pro:
fail += 1
self.__fail_num += fail
self.__cure_num += number - fail
self.__expect_pro = self.__cure_num / self.__patient_num
return (number-fail, fail)
def sampling_pro(self):
"""通过汤普森采样，获得guess_pro"""
random_variable = scipy.stats.beta(self.cure_num,
self.fail_num)
self.guess_pro = random_variable.rvs() #采样
class Pills:
"""
>>> pros = [i / 10 for i in range(1, 9)]
>>> pills = Pills(pros, 2, 2)
>>> pills.pills_num
8
"""
def __init__(self, pros, cure_num, fail_num):
"""
初始化
:param pros:每个药丸的治疗概率
:param cure_num: 先验a
:param fail_num: 先验b
"""
self.__pills = []
for pill_pro in pros:
pill = Pill(pill_pro, cure_num, fail_num)
self.__pills.append(pill)
self.guess_pros = [pill.guess_pro
for pill in self.__pills]
self.__pills_num = len(self.__pills)
self.__patient_nums = [pill.patient_num
for pill in self.__pills]
self.__expect_pros = [pill.expect_pro
for pill in self.__pills]
self.process = []
@property
def pills_num(self):
return self.__pills_num
@property
def expect_pros(self):
"""获取期望概率"""
return self.__expect_pros
@property
def patient_nums(self):
return self.__patient_nums
def sampling_pros(self):
"""为药丸们采样"""
for pill in self.__pills:
pill.sampling_pro()
def update_pros(self):
"""更新概率"""
self.guess_pros = [
pill.guess_pro
for pill in self.__pills
]
self.__expect_pros = [
pill.expect_pro
for pill in self.__pills
]
def choose_pill(self):
"""选择药丸
实际上，就是寻找治疗概率最大的那个
"""
self.sampling_pros()
self.update_pros()
pros = np.array(self.guess_pros)
location = pros.argmax()
pill = self.__pills[location]
return (pill, location)
def curing(self, number=1):
"""选择药丸，进行治疗,只针对一种药丸"""
pill, location = self.choose_pill()
cure, fail = pill.curing(number)
self.process = self.process + [location] * number #记录下过程
self.patient_nums[location] = pill.patient_num
return cure, fail
def all_curing(self, number=1):
"""所有药丸都要试一遍"""
self.update_pros()
for k, pill in enumerate(self.__pills):
pill.curing(number)
self.patient_nums[k] = pill.patient_num
import matplotlib.pyplot as plt
pros = [i / 10 for i in range(1, 9)]
pills = Pills(pros, 2, 2)
opacity = 0.4
bar_width = 0.08
for i in range(100):
pills.all_curing(3)
expect_pros = pills.expect_pros
nums = pills.patient_nums
fig, ax = plt.subplots()
ax.bar(pros, expect_pros, bar_width,
alpha=opacity, color='r')
ax.set_title("Beta-Bernoulli")
ax.set_ylabel("Expected probability")
for j in range(8):
ax.text(pros[j] - 0.02, expect_pros[j] - 0.03, round(expect_pros[j], 2))
fig.savefig("ben/pic{0}".format(
i
))
plt.close()
print(pills.guess_pros)
print(pills.patient_nums)
if __name__ == "__main__":
import doctest
doctest.testmod()
# 高斯过程
import numpy as np
PIC = 0
class BayesOpti_GP1:
"""
贝叶斯优化，基于高斯过程
只是用于一维的，可以进行扩展
"""
def __init__(self, x, y, theta, sigma):
"""
:param x: 位置坐标
:param y: 输出
:param theta: 包括theta0 和 theta1
:param sigma:
"""
self.__x = [x] # x是ndarray
self.y = [y] #y也是ndarray
self.theta0 = theta[0]
self.theta1 = theta[1] #尺度参数
self.sigma = sigma[0]
self.n = len(self.__x) #即x的个数
self.I = np.diag(np.ones(self.n, dtype=float)) #单位矩阵
self.minimum = min(self.y) #最小值 用于EI
@property
def x(self):
"""获取属性x"""
return self.__x
"""
我们设定x为私有变量，所以原则上允许改变x
def set_x(self, x):
self.__x.append(x)
"""
def add_y(self, y):
"""添加y同时更新最小值"""
self.y.append(y)
self.minimum = min(self.y)
def kernel(self, r2):
"""为了简单,采用exp kernel"""
return self.theta0 ** 2 * np.exp(-1/2 * r2)
def matrix_K(self):
"""更新矩阵K，同时还有一系列衍生的矩阵，
实际上有很多计算的浪费在此处，但是不想纠结这么多了
"""
self.K = np.zeros((self.n, self.n), dtype=float)
self.m_grad1 = np.zeros((self.n, self.n), dtype=float)
for i in range(self.n):
for j in range(i, self.n):
r2 = (self.__x[i]-self.__x[j]) ** 2
* self.theta1 ** 2
self.K[i, j] = self.kernel(r2)
self.m_grad1[i, j] = -self.K[i, j] * r2 / self.theta1 #关于theta1的导数
self.K[j, i] = self.K[i, j]
self.m_grad1[j, i] = self.m_grad1[i, j]
self.K_sigma = self.K + self.sigma ** 2 * self.I#K+sigma^2I
self.K_sigma_inverse = np.linalg.inv(self.K_sigma) # 逆
self.K_sigma_det = np.linalg.det(self.K_sigma) #行列式
self.m_grad0 = 2 * self.K / self.theta0 #关于theta0的导数
def grad_matrix(self, label):
"""根据需要选择导数矩阵。。。"""
if label is 0:
matrix = self.m_grad0
elif label is 1:
matrix = self.m_grad1
else:
matrix = self.I * 2 * self.sigma
return matrix
def grad_detA(self, label):
"""对行列式求导"""
matrix = self.grad_matrix(label)
grad = 0.
for i in range(self.n):
ksigma = self.K_sigma.copy()
ksigma[i] = matrix[i]
grad += np.linalg.det(ksigma)
y = np.array(self.y, dtype=float)
grad = (1. - y @ self.K_sigma_inverse @ y) * grad
return grad
def grad_A(self, label):
"""part1求导，不知道该怎么描述"""
A_star = np.zeros((self.n, self.n), dtype=float)
def grad_A_ij(self, matrix, i, j):
m1 = np.delete(matrix, j, axis=0)
m1 = np.delete(m1, i, axis=1)
m2 = np.delete(self.K_sigma, j, axis=0)
m2 = np.delete(m2, i, axis=1)
grad = 0.
for k in range(self.n-1):
m3 = m2.copy()
m3[k] = m1[k]
grad += np.linalg.det(m3)
return grad * (-1) ** (i + j)
matrix = self.grad_matrix(label)
for i in range(self.n):
for j in range(i, self.n):
A_star[i, j] = grad_A_ij(self, matrix, i, j)
A_star[j, i] = A_star[i, j]
y = np.array(self.y, dtype=float)
return y @ A_star @ y
def grad_pra(self, label):
"""落实道每个参数的导数，这回是真的
导数了，而不是中间的过渡的矩阵
"""
part1 = self.grad_A(label)
part2 = self.grad_detA(label)
grad = (part1 + part2) / self.K_sigma_det
return grad
def marginal_y(self):
"""y的边际分布，省略了很多东西，负号也去了，
因为用的是梯度下降"""
y = np.array(self.y, dtype=float)
return y @ self.K_sigma_inverse @ y
+ np.log(self.K_sigma_det)
def update_pra(self):
"""更新参数，如果n为1是一种特殊情况
采用梯度下降方法，而且是很愚蠢的那种
"""
if self.n is 1:
self.matrix_K()
self.theta0 = self.y[0] * np.sqrt(2) / 2
self.sigma = self.y[0] * np.sqrt(2) / 2
return 1
grad = [999., 999., 999.]
t = 1
self.matrix_K()
min = self.marginal_y()
min_pra = [self.theta0, self.theta1, self.sigma]
while True:
if sum(list(map(abs, grad))) < 1e-4 or t > 200:
break
grad[0] = self.grad_pra(0)
grad[1] = self.grad_pra(1)
grad[2] = self.grad_pra(2)
step = max([0.013, 1/t]) #学习率
self.theta0 = self.theta0 - step * grad[0]
self.theta1 = self.theta1 - step * grad[1]
self.sigma = self.sigma - step * grad[2]
self.matrix_K()
y = self.marginal_y()
if y < min:
min = y
min_pra = [self.theta0, self.theta1, self.sigma]
else:
pass
t += 1
self.theta0 = min_pra[0]
self.theta1 = min_pra[1]
self.sigma = min_pra[2]
return 1
def find_newx(self):
"""根据PI寻找x"""
x = np.linspace(0, 1, 1000)
y = np.array(self.y)
z = []
u = []
for item in x:
k = []
for xi in self.__x:
r2 = (xi - item) ** 2 * self.theta1 ** 2
k.append(self.kernel(r2))
k = np.array(k)
q = (self.minimum - k @ self.K_sigma_inverse @ y) /
(self.theta0 ** 2 - k @ self.K_sigma_inverse @ k)
u.append(k @ self.K_sigma_inverse @ y)
z.append(q)
z = np.array(z)
u = np.array(u)
import matplotlib.pyplot as plt
#plt.cla() #清空当前axes
#plt.plot(x, u)
#plt.pause(0.1)
plt.cla()
fig, ax = plt.subplots()
ax.plot(x, u)
fig.savefig("bayes/pic{0}".format(
PIC
))
newx = x[z.argmax()]
self.__x.append(newx)
self.n += 1
self.I = np.diag(np.ones(self.n, dtype=float))
return newx
def f(x): #实际的函数
return (x - 0.3) ** 2 + 0.2 * np.sin(20 * x)
def f2(x): #加了噪声的函数
return (x - 0.3) ** 2 + 0.2 * np.sin(20 * x)
+ np.random.randn() * 0.05
points = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
count = 3
x0 = points[count]
y0 = f(x0)
theta = np.random.randn(2) * 100 #给定初始的参数 随机给的
sigma = np.random.randn(1) * 100
test = BayesOpti_GP1(x0, y0, theta, sigma)
for i in range(10): #进行10次评估
test.update_pra()
newx = test.find_newx()
newy = f2(newx)
test.add_y(newy)
PIC += 1
x = test.x
y = test.y
x1 = np.linspace(0, 1, 1000)
y1 = [f(item) for item in x1]
import matplotlib.pyplot as plt
plt.cla()
print(x, y)
print(test.theta0, test.theta1, test.sigma)
fig, ax = plt.subplots()
ax.plot(x1, y1)
ax.scatter(x, y)
for i in range(len(x)):
plt.text(x[i], y[i]+0.02, str(i+1), size=10)
plt.show()
附录
Beta-Bernoulli 模型的推导
又(p(D|mathrm{w}_a)=w_a^{n_{a,1}}(1-w_a)^{n_{a,0}})
所以，(p(D|mathrm{w}_a)p(mathrm{w}_a))的核为(w_a^{n_{a,1}+alpha-1}(1-w_a)^{n_{a,0}+eta-1}),满足的是(mathrm{Beta}(w_a|a+n_{a,1},eta+n_{a,0}))
证毕。
线性模型 后验分布(p(mathrm{w},sigma^2|D_n))的推导
只需考虑(p(D_n|mathrm{w},sigma^2)p(mathrm{w},sigma^2))的核即可。
先配(mathrm{w})，再配第二部分，即可得：
(y)边际分布推导
如果我们令(f=mathrm{Xw})，那么(p(y|mathrm{X},sigma^2))也可以表示为：
这个积分比先前的积分要容易求解（至少前面的那个我直接求不出来）。
依旧采用核方法：
其中:
上面第二个(propto)的前半部分在积分中相当于常数可以提出，后半部分会被积分掉。所以，现在我们只需要证明：
又
容易证明((A+B)^{-1}=A^{-1}(A^{-1}+B^{-1})^{-1}B^{-1})（(A,B)可逆）
所以，
令(C = mathrm{I}-Sigma(Sigma+sigma^2)^{-1}),
等式俩边，右乘((Sigma+sigma^2))得：
所以，
所以，
得证，
证毕.
(f_*)的后验分布
我们不加推导地给出：
根据高斯过程的性质，存在如下的联合分布：
其中，(mathrm{X}_*)表示预测输入，而(f_*)表示预测输出，(m_*=u_0(mathrm{X_*}))，(K_*^mathrm{T}={k(mathrm{x}_1,mathrm{X_*}), ldots, k(mathrm{x}_n,mathrm{X}_*)}),(K_{**}=k(mathrm{X}_*, mathrm{X}_*))
我们有：
所以，同样的，我们只需要考虑分子关于$f_* $的核就行了。
其中:
容易推得(f_*)的均值(mu_*)和协方差(sigma_*^2)为：
根据Schur补和分块矩阵求逆（凸优化-p622）的性质，我们可以得到：
所以，我们得到：
证毕.