贝叶斯优化

[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.

简介

贝叶斯优化-marsggbo

首先，贝叶斯优化能干什么？给我的感觉是无所不能，当然其效果有些可能不尽如人意。贝叶斯优化，可以做回归的东西（虽然总感觉这个东西只是一个附属品），然而主要是去解决一个“优化问题”。
贝叶斯优化解决的是下面类型的问题：
$${\mathrm{x}}^{\ast }=\underset{\mathrm{x}\in \mathcal{X}}{\mathrm{argmax}}f(\mathrm{x})$$
注: 使用"argmin"并无实质上的不同，事实上[1]中采用的便是"argmin"。
往往，$f$我们并不知道，所以，这类问题很难采用经典的梯度上升("argmin"则梯度下降)来解决。贝叶斯优化采用概率代理模型来应对。$\mathrm{x}$是决策，往往称$\mathcal{X}$为决策空间。药物配方是一种决策，神经网络卷积核大小等也可以看成一种决策。而且，这种决策与最后的输出的关系（即$f$）往往很难知晓。这也正是贝叶斯优化的强大之处。

贝叶斯优化框架

上面俩幅图分别来自[2]和[1]，因为一些符号的差异，往下除特别指明，采用的均为[2]中的符号。
贝叶斯优化，每一次迭代，首先在代理模型的“先验”下，通过最大化采集函数（该函数往往是对评估点的分布以及$f(x)$的提升的一种权衡(trade-off)）。新的评估点，作为输入传入系统，获得新的输出，以此来更新$D$和概率代理模型。
其中 D t = { ( x 1 , y 1 ) , … , ( x t , y t ) } D_t=\{(\mathrm{x_1}, y_1), \ldots, (\mathrm{x_t},y_t)\} Dt​={(x1​,y1​),…,(xt​,yt​)}

上面这幅图，便是贝叶斯优化的一个简单演示。黑色虚线表示目标函数$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(\mathrm{w}|D)=\frac{p(D|\mathrm{w})p(\mathrm{w})}{p(D)}$$
现在问题来呢，我们还不知道$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\dots K$来表示药物，$w_a$表示第$a$种药物成功治愈病患的可能性，而 y i ∈ { 0 , 1 } y_i \in \{0, 1\} yi​∈{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{B}\mathrm{e}\mathrm{t}\mathrm{a}$分布作为$\mathrm{w}$的先验分布，因为这是其共轭先验分布。
$$p(\mathrm{w}|\alpha ,\beta )=\prod _{a=1}^K\mathrm{B}\mathrm{e}\mathrm{t}\mathrm{a}(w_a|\alpha ,\beta )$$
定义：
$$\begin{array}{l}{n}_{a,0}=\sum _{i=1}^n\mathbb{I}(y_i=0,a_i=a)\\ n_{a,1}=\sum _{i=1}^n\mathbb{I}(y_i=1,a_i=a)\end{array}$$
其中$n_{a,0}$表示$n$次评估中，选中$a$药物且治疗失败的数目，$n_{a,1}$则反之。$\mathbb{I}(\cdot )$只有$(\cdot )$成立为1否则为0。
那么，$\mathrm{w}$的后验概率为：
$$p(\mathrm{w}|D)=\prod _{a=1}^K\mathrm{B}\mathrm{e}\mathrm{t}\mathrm{a}(w_a|\alpha +n_{a,1},\beta +n_{a,0})$$
上述推导见附录。
从上述也能发现，超参数$\alpha ,\beta$表示的治愈数和失败数。下图是以$\mathrm{B}\mathrm{e}\mathrm{t}\mathrm{a}(w|2,2)$为先验的一个例子。

汤普森采样-wiki
那么在$D_n$的基础上，如果找$a_{n+1}$呢。以下采用的是汤普森采样（或后验采样）：
$${\alpha }_{n+1}=\underset{a}{\mathrm{argmax}}{f}_{\tilde{\mathrm{w}}}(a)$$
$\tilde{\mathrm{w}}\sim p(\mathrm{w}|D_n)$,即$\tilde{\mathrm{w}}$从$\mathrm{w}$的后验分布中采样得到。
该模型的好处是：

• 只有$\alpha ,\beta$俩个超参数
• 是对exploration和exploitation的一种比较好的权衡
• 比较容易和蒙特卡洛采样-HFUT_qianyang结合
• 汤普森采样的随机性使其容易推广到批量处理（不懂）

下面是该模型的算法：

线性模型(Linear models)

上述的模型在应对组合类型的时候会显得捉襟见肘。比方，我们在从$d$个元素中寻求一种搭配，每个元素有 { 0 , 1 } \{0, 1\} {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}$成了权重向量。这只是代理模型的一部分，因为并没有体现“概率”的部分。
$$y_a\sim N({\mathrm{x}}_a^{\mathrm{T}}\mathrm{w},{\sigma }^2)$$
组合$a$的观测值为$y_a$，服从正态分布。很自然地，我们同样选择共轭先验分布作为$\mathrm{w},{\sigma }^2$的先验分布：normal_inverse_gamma-wiki
N I G ( w , σ 2 ∣ w 0 , V 0 , α 0 , β 0 ) = ∣ 2 π σ 2 V 0 ∣ − 1 2 e x p { − 1 2 σ 2 ( w − w 0 ) T V 0 − 1 ( w − w 0 ) } × β 0 α 0 Γ ( α 0 ) ( σ 2 ) α + 1 e x p { − β 0 σ 2 } \begin{array}{l} \mathrm{NIG}(\mathrm{w}, \sigma^2|\mathrm{w}_0,V_0,\alpha_0,\beta_0)=\\ |2\pi \sigma^2 V_0|^{-\frac{1}{2}}\mathrm{exp}\{-\frac{1}{2\sigma^2}(\mathrm{w-w_0})^\mathrm{T}V_0^{-1}(\mathrm{w-w_0})\} \times \frac{\beta_0^{\alpha_0}}{\Gamma(\alpha_0)(\sigma^2)^{\alpha + 1}}\mathrm{exp}\{-\frac{\beta_0}{\sigma^2}\} \end{array} NIG(w,σ2∣w0​,V0​,α0​,β0​)=∣2πσ2V0​∣−21​exp{−2σ21​(w−w0​)TV0−1​(w−w0​)}×Γ(α0​)(σ2)α+1β0α0​​​exp{−σ2β0​​}​
$\mathrm{N}\mathrm{I}\mathrm{G}$分布有4个超参数，而$\mathrm{w},{\sigma }^2$的后验分布（$D_n$的条件下）满足：

其中$\mathrm{X}$的第$i$行为${\mathrm{x}}_{{\alpha }_i}$。
推导见附录。

关于${\alpha }_{n+1}$的选择，同样可以采用汤普森采样：
$${\alpha }_{n+1}=\underset{a}{\mathrm{argmax}}{\mathrm{x}}_a^{\mathrm{T}}\tilde{\mathrm{w}}$$
其中$\tilde{\mathrm{w}}\sim p(\mathrm{w}|D_n)$

线性模型有很多扩展：
$f(x)=\Phi (\mathrm{x}{)}^{\mathrm{T}}\mathrm{w}$
$f(x)=g({\mathrm{x}}^{\mathrm{T}}\mathrm{w})$

其中，$\Phi (\mathrm{x})=({\varphi }_1(\mathrm{x}),\dots ,{\varphi }_J(\mathrm{x}){)}^{\mathrm{T}}$，${\varphi }_j(\mathrm{x})$常常为：
e x p { − ( x − z j ) T Λ ( x − z j ) 2 } \mathrm{exp}\{ -\frac{(\mathrm{x-z_j})^{\mathrm{T}}\Lambda(\mathrm{x-z_j})}{2} \} exp{−2(x−zj​)TΛ(x−zj​)​}
和
e x p { − i x T w j } \mathrm{exp}\{ -i\mathrm{x^T}w_j\} exp{−ixTwj​}
这里，$\Lambda$， { z j } j = 1 J \{\mathrm{z}_j \}_{j=1}^{J} {zj​}j=1J​， { w j } j = 1 J \{ w_j\}_{j=1}^{J} {wj​}j=1J​均为超参数，至于这些超参数怎么更新，我不大清楚。

非参数模型

非参数模型不是指没有参数，而是指参数（数量）不定。

我们先来看如何把先前的线性模型转换成非参数模型。
我们假设${\sigma }^2$是固定的，且$p(\mathrm{w}|V_0)=\mathcal{N}(0,V_0)$,即服从均值为$0$,协方差矩阵为$V_0$的多维正态分布。那么，$y\sim \mathcal{N}(\mathrm{X}\mathrm{w},{\sigma }^2\mathrm{I})$，我们可以积分掉$\mathrm{w}$从而获得$y$的一个边际分布:

推导见附录。
就像先前已经提到过的，我们可以引入$\Phi =({\varphi }_i,\dots ,{\varphi }_J{)}^T$,
将资料（设计）矩阵$X$映射到${\mathbb{R}}^{n\times J}$，如此一来，相应的边际分布也需发生变化：

注意到$\Phi V_0{\Phi }^T$，事实上，我们不需要特别指明$\Phi$，而只需通过kernel.

$K$将是$\Phi V_0{\Phi }^T$的一个替代。采用这个策略，比原先在计算上和可解释性上更有优势。
不过，还有另外一个问题，如果去寻找下一个评估点呢。寻找下一个评估点，需要我们做预测，但是上面的边际分布显然是无法进行预测的。不过，我们可以通过条件公式来获得：

${\mathrm{X}}_{\ast }$是新的位置,而$y_{\ast }$是相应的预测，二者都可以是向量。
分子部分是一个联合的高斯分布。到此，我们实际上完成了一个简易的高斯过程，下面正式介绍高斯过程。

高斯过程

高斯过程-wiki
高斯过程-火星十一郎

高斯过程$f(\mathrm{x})\sim GP({\mu }_0,k)$，其核心便是均值函数${\mu }_0$(在贝叶斯优化中，我们常常选择其为0)和协方差函数$k({\mathrm{x}}_i,{\mathrm{x}}_j)$，而观测值$y=f+\epsilon$。通过高斯过程得到的序列$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$是随机变量$\epsilon$的方差。
于是，我们可以像之前所做的那样，求边际分布，和$y_{\ast }$的分布。
首先，
$$p(y|\mathrm{X},{\sigma }^2)=\mathcal{N}(y|m,K+{\sigma }^2\mathrm{I})$$
我们并没有给出这个证明，因为这个不难验证。接下来，为了预测，我们需要求后验分布。论文此时并没有选择$y_{\ast }$，而是$f_{\ast }$的后验分布，这点倒是可以理解，比较我们的目标就是最大化$f(\mathrm{x})$。不过，论文给出的是标量（也就是只有一个预测值），实际上，很容易扩展到多个，在附录里，我们给出多个的后验分布的推导。
我们先给出一个的后验分布，依旧是正态分布：

常用的一些kernels

Kernel method - wiki
Matern covariance function - wiki

文献[1]给出的形式如下(实际上是$d=1$的情况)：

其中，$r=|x-x^{\prime }|$，$v$为平滑参数，$l$为尺度参数，$K_v$为第二类变形贝塞尔函数。
同时给出了几种常用的Matern协方差函数。

文献[2]给出的是另外一种表示方式：

其中，$r^2=(\mathrm{x}-{\mathrm{x}}^{\prime }{)}^{\mathrm{T}}\mathrm{\Lambda }(\mathrm{x}-{\mathrm{x}}^{\prime })$，$\Lambda$是一个对角矩阵，其对角线元素为${\theta }_i^2,i=1,\dots ,d$。

这些参数可以这么理解：

• $v$，被称为平滑参数，因为，从使用Matern协方差函数的高斯过程中采样的目标函数$f(x)$是$\lfloor v-1\rfloor$次均方可微的（为什么？）。
• $l$和$\theta$的功能相似，反映该维度的重要性（前者是越小越重要，后者越大越重要）。

上面的一些参数，会在下面给出一些更新的方法。

边际似然

log 边际似然函数可以表示为：

从图中可以看到，等式右边被分成了三部分，三者有不同含义：

• 第一项，表示模型和数据的拟合程度的好坏，以马氏距离为指标；
• 第二项，表示模型的复杂度；
• 第三项，是数据点$n$的线性函数，表示数据的不确定性随着数据增大而减少（？）。

一个非常自然的想法是，对上述似然函数进行极大似然估计，从而获得$\theta :=({\theta }_{0:d},{\mu }_0,{\sigma }^2)$的估计。

复杂度

每一次高斯过程的复杂度在$O(n^3)$级别左右，这是由计算矩阵的逆所带来的。通过Cholesky分解，可以降为$O(n^2)$。
所以产生了一些算法，试图克服这个困难。

sparse pseudo-input Gaussian processes (SPGP)

SPGP从n个输入中选择m个伪输入替代，从而达到降秩的目的。同时，这m个伪输入的位置也作为参数（虽然我不知道怎么去更新）。好处自然是，
能够把复杂度降为$0(n{m}^2+m^3)$。缺点是，参数相对比较多，容易产生“过拟合”现象。

Sparse spectrum Gaussian processes(SSGP)

由Bochner定理得，任何稳定得kernel都有一个正定得傅里叶谱表示：

之后，通过蒙特卡洛方法，采样m个样本频率，来近似估计上诉的积分。从而获得近似的协方差函数，当数据集较小时，SSGP同样易产生“过拟合”现象。

随机森林

随机森林 - Poll的笔记

随机森林可以作为高斯过程的一种替代。缺点是，数据缺少的地方，估计的并不准确（边际更是常数）。另外，由于随机森林不连续，也就不可微，所以无法采用梯度下降（或上升）的方法来更新参数。另外不解的是，随机森林的参数，即便我们给了一个先验分布，可其后验分布如何求呢？

采集函数

首先我们有一个效用函数$U:{\mathbb{R}}^d\times \mathbb{R}\times \Theta \to \mathbb{R}$，顾名思义，效用函数，是反映评估$\mathrm{x}$和对应的函数值$v=f(\mathrm{x})$（在$\theta$条件下）的一个指标。论文[1]并没有引入这个效用函数，论文[2]引入这个概念应该是为了更好的说明。

一种采集函数的选择，便是期望效用：

其实蛮奇怪的，因为对$v$求期望也就罢了，这个采集函数对$\theta$也求了期望，我的理解是，这样子更加“模糊”了，如果选择极大似然等方式产出的“精准”的$\theta$，或许不能够很好的让评估点足够分散，而陷入局部最优。而且，这样子做，似乎就没有必要去估计参数$\theta$，虽然代价是求期望。

从下面的一些算法中我们可以看到，往往没有${\mathbb{E}}_{\theta }(\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})$，令其为$\tau$。往往，$\tau ={\min }_{\mathrm{x}\in \mathrm{X}}f(\mathrm{x})$。
采集函数为：

注意，这里的$\Phi$是标准正态分布的概率函数。
这个采集函数里的效用函数是：
$$U(\mathrm{x},v,\theta )=\mathbb{I}[v>\tau ]$$
其中$\mathbb{I}$为指示函数。
当$\tau$就是$f(\mathrm{x})$的最小值时，PI的效果非常好。
PI一个“弊端”是，只在乎提升的概率，而在乎提升的幅度，而，EI就涵盖了这俩方面。

EI(expected improvement)

通常，其提升函数由下式表示：

而相应的的采集函数是：

其中$\varphi$是标准正态分布的概率密度函数。式子通过积分变量替换可以推得。
实际上$I$就是效用函数$U$。

UCB(upper confidence bound)

采集函数为：

这个采集函数，可以这么理解，对于任意一个$\mathrm{x}$，它有一个均值${\mu }_n(\mathrm{x})$，有一个标准差${\sigma }_n(\mathrm{x})$（体现浮动范围和程度），${\beta }_n{\sigma }_n(\mathrm{x})$我们认为比较可靠的界，我们认为，$f(\mathrm{x})$有较大可能达到${\mu }_n(\mathrm{x})+{\beta }_n{\sigma }_n(\mathrm{x})$的值。所以最大化采集函数，就是最大化我们的这一种希望。
论文[2]中说$\beta$的选择往往是Chernoff-Hoeffding界。听起来很玄乎，但是，UCB现在貌似非常火。另外，有一套理论，能够引导和规划超参数${\beta }_n$，使得能够达到最优。

基于信息的策略

不同之前的策略，基于信息的策略，依赖全局最优解${\mathrm{x}}^{\ast }$的后验分布$p_{\ast }(\mathrm{x}|D_n)$。该分布，隐含在函数$f$的后验分布里（不同的$f$有不同的全局最优解，从而${\mathrm{x}}^{\ast }$也有一个后验分布）。

熵搜索算法旨在寻找能够极大减少不确定度的评估点${\mathrm{x}}^{\ast }$，从另外一个角度来说，${\mathrm{x}}^{\ast }$给与了我们最多的信息。

因此效用函数为：

此时的采集函数为：

其中$H({\mathrm{x}}^{\ast }|{\mathrm{D}}_{\mathrm{n}})=\int p({\mathrm{x}}^{\ast }|{\mathrm{D}}_{\mathrm{n}})\log p({\mathrm{x}}^{\ast }|D_n)d{\mathrm{x}}^{\ast }$代表微分熵，$y\sim \mathcal{N}({\mu }_n(\mathrm{x}),{\sigma }_n^2(\mathrm{x})+{\sigma }^2)$。
不过，估计$p(\mathrm{x}\ast |D_n)$或者$H$都绝非易事，计算量也让人望而却步。

一种新的方法PES(predictive entropy search)很好地改善了这些问题。
利用互信息的对称性（不懂），我们可以将${\alpha }_{ES}(\mathrm{x})$改写为：

使得我们不必苦于$p({\mathrm{x}}^{\ast }|D_n)$上。

采集函数的组合

不同采集函数的组合或许能够达到更好的特性。
一种方案是，不同采集函数会给出一系列最优评估点的候选。通过一个准则来判断孰优孰劣。
ESP（entropy search portfolio）采用的准则如下：

即选择那个让不确定性最小的候选点。

出于实际的考量

处理超参数

处理（优化）超参数一般有如下的策略：

• 通过$y$的边际分布，进行极大似然估计，获得${\hat{\theta }}_n^{ML}$
• 给参数$\theta$一个先验分布$p(\theta )$，通过贝叶斯公式获得其后验分布$p(\theta |D_n)$，可得最大后验估计${\hat{\theta }}_n^{MAP}$
• 同样用到后验分布$p(\theta |D_n)$，对$\theta$进行边际化处理：
  $${\alpha }_n(\mathrm{x})={\mathbb{E}}_{\theta }\alpha (\mathrm{x};\theta )=\int p(\theta |D_n)\alpha (\mathrm{x};\theta )\mathrm{d}\theta$$ 这部分不必再估计参数$\theta$，代价是需要估计一个期望，这项工作有不同的方案：蒙特卡洛技术，贝叶斯蒙特卡洛技术等。

采集函数的优化

实际上，采集函数的优化并不简单，因为其往往非凸。有基于梯度的方法（如果可微的话），还有网格搜索，自适应协方差矩阵进化策略，多启动局部搜索等。
最近还流行用空间切分的方法替代贝叶斯模型，如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个无差别病患。下图是一个例子：

再下面的例子是，每次试验不选择药丸，统统进行试验：

高斯过程

选择最简单的一维高斯过程，超参数有${\theta }_0,{\theta }_1,\sigma$,选择的kernel为 θ 0 2 e x p { − 1 2 r 2 } \theta_0^2 exp\{-\frac{1}{2}r^2 \} θ02​exp{−21​r2}，采集函数为$EI$。优化超参数使用梯度下降（因为别的方法不怎么会），优化采集函数使用的是网格搜索。另外，我们给输出附加了方差为0.0001（基本上没有）的白噪声。

采用[1]中的例子：

首先，是在我们不知道$\tau$的情况下，我们取初始值0.1，0.2， 0.4，0.5，0.7，0.9.，上面的图是选取11个点后的均值函数，下面的是真实的曲线，上面的点是每次选取的点。

接下来，我们固定$\tau =-0.2$

最后再来一个加入方差为0.0025白噪声的$\tau =-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 模型的推导

$$\mathrm{B}\mathrm{e}\mathrm{t}\mathrm{a}(w_a|\alpha ,\beta )=\frac{\Gamma (\alpha +\beta )}{\Gamma (\alpha )\Gamma (\beta )}{w}_a^{\alpha -1}(1-w_a{)}^{\beta -1},w_a\in [0,1]$$
$$\begin{array}{ll}p(\mathrm{w}|D)& =\prod _{a=1}^{k}p({\mathrm{w}}_a|D)\\ & \propto \prod _{a=1}^{k}p(D|{\mathrm{w}}_a)p({\mathrm{w}}_a)\end{array}$$
又$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}+\beta -1}$,满足的是$\mathrm{B}\mathrm{e}\mathrm{t}\mathrm{a}(w_a|a+n_{a,1},\beta +n_{a,0})$
证毕。

线性模型 后验分布$p(\mathrm{w},{\sigma }^2|D_n)$的推导

只需考虑$p(D_n|\mathrm{w},{\sigma }^2)p(\mathrm{w},{\sigma }^2)$的核即可。

p ( D n ∣ w , σ 2 ) ∝ 1 σ n e x p { − ( y − X w ) T ( y − X w ) 2 σ 2 } p(D_n|\mathrm{w},\sigma^2)\propto \frac{1}{\sigma^n}\mathrm{exp}\{-\frac{(y-\mathrm{Xw})^{\mathrm{T}}(y-\mathrm{Xw})}{2\sigma^2}\} p(Dn​∣w,σ2)∝σn1​exp{−2σ2(y−Xw)T(y−Xw)​}
p ( w , σ 2 ) ∝ 1 σ 2 α 0 + 3 e x p { − ( w − w 0 ) T V 0 − 1 ( w − w 0 ) + 2 β 0 2 σ 2 } p(\mathrm{w}, \sigma^2)\propto \frac{1}{\sigma^{2\alpha_0+3}}\mathrm{exp}\{-\frac{(\mathrm{w-w_0})^{\mathrm{T}}V_0^{-1}(\mathrm{w-w_0})+2\beta_0}{2\sigma^2} \} p(w,σ2)∝σ2α0​+31​exp{−2σ2(w−w0​)TV0−1​(w−w0​)+2β0​​}
先配$\mathrm{w}$，再配第二部分，即可得：

$y$边际分布推导

如果我们令$f=\mathrm{X}\mathrm{w}$，那么$p(y|\mathrm{X},{\sigma }^2)$也可以表示为：
$$\begin{array}{ll}p(y|\mathrm{X},{\sigma }^2)& =\int p(y|f,{\sigma }^2)p(f|\mathrm{X},{\mathrm{V}}_0)df\\ & =\int \mathcal{N}(y|f,{\sigma }^{2}I)\mathcal{N}(f|0,\mathrm{X}{\mathrm{V}}_0{\mathrm{X}}^{\mathrm{T}})df\end{array}$$
这个积分比先前的积分要容易求解（至少前面的那个我直接求不出来）。
依旧采用核方法：
p ( y ∣ X , σ 2 ) ∝ e x p { − ( y − f ) T ( σ 2 I ) − 1 ( y − f ) + f T ( X V 0 X T ) − 1 f 2 } ∝ e x p { − y T ( ( σ − 2 − σ − 2 ( Σ − 1 + σ − 2 ) − 1 σ − 2 ) y 2 } × e x p { − ( f − μ ) T ( Σ − 1 + σ − 2 ) ( f − μ ) 2 } \begin{array}{ll} p(y|\mathrm{X}, \sigma^2) & \propto \mathrm{exp} \{-\frac{(y-f)^{\mathrm{T}}(\sigma^2 I)^{-1}(y-f)+f^{\mathrm{T}}(\mathrm{XV_0X^T})^{-1}f}{2}\}\\ & \propto \mathrm{exp} \{-\frac{y^{\mathrm{T}}((\sigma^{-2}-\sigma^{-2}(\Sigma^{-1}+\sigma^{-2})^{-1}\sigma^{-2})y}{2}\} \\ & \quad \times \mathrm{exp}\{ -\frac{(f-\mu)^{\mathrm{T}}(\Sigma^{-1}+\sigma^{-2})(f-\mu)}{2}\} \end{array} p(y∣X,σ2)​∝exp{−2(y−f)T(σ2I)−1(y−f)+fT(XV0​XT)−1f​}∝exp{−2yT((σ−2−σ−2(Σ−1+σ−2)−1σ−2)y​}×exp{−2(f−μ)T(Σ−1+σ−2)(f−μ)​}​
其中:
$$\begin{array}{l}\mu =({\Sigma }^{-1}+{\sigma }^{-2}{)}^{-1}{\sigma }^{-2}y\\ \Sigma =\mathrm{X}{\mathrm{V}}_0{\mathrm{X}}^{\mathrm{T}}\end{array}$$
上面第二个$\propto$的前半部分在积分中相当于常数可以提出，后半部分会被积分掉。所以，现在我们只需要证明：
$$(({\sigma }^{-2}-{\sigma }^{-2}({\Sigma }^{-1}+{\sigma }^{-2}{)}^{-1}{\sigma }^{-2}{)}^{-1}=\Sigma +{\sigma }^2\mathrm{I}$$
又
$$\begin{array}{ll}(({\sigma }^{-2}-{\sigma }^{-2}({\Sigma }^{-1}+{\sigma }^{-2}{)}^{-1}{\sigma }^{-2})& =\frac{{\sigma }^2-({\Sigma }^{-1}+{\sigma }^{-2}{)}^{-1}}{{\sigma }^2}\end{array}$$
容易证明$(A+B{)}^{-1}=A^{-1}(A^{-1}+B^{-1}{)}^{-1}{B}^{-1}$（$A,B$可逆）
所以，
$$({\Sigma }^{-1}+{\sigma }^{-2}{)}^{-1}=\Sigma (\Sigma +{\sigma }^2{)}^{-1}{\sigma }^2$$
令$C=\mathrm{I}-\Sigma (\Sigma +{\sigma }^2{)}^{-1}$,
等式俩边，右乘$(\Sigma +{\sigma }^2)$得：
$$C(\Sigma +{\sigma }^2)={\sigma }^2$$
所以，
$$C=\mathrm{I}-\Sigma (\Sigma +{\sigma }^2{)}^{-1}={\sigma }^2(\Sigma +{\sigma }^2{)}^{-1}$$
所以，
$$(({\sigma }^{-2}-{\sigma }^{-2}({\Sigma }^{-1}+{\sigma }^{-2}{)}^{-1}{\sigma }^{-2})=(\Sigma +{\sigma }^2\mathrm{I}{)}^{-1}$$
得证，
证毕.

$f_{\ast }$的后验分布

我们不加推导地给出：
$$p(y|\mathrm{X},{\sigma }^2)=\mathcal{N}(y|m,K+{\sigma }^2\mathrm{I})$$
根据高斯过程的性质，存在如下的联合分布：
$$\left[\begin{array}{c}y\\ f_{\ast }\end{array}\right]\sim \mathcal{N}(\left[\begin{array}{c}m\\ m_{\ast }\end{array}\right],\left[\begin{array}{cc}K+{\sigma }^2\mathrm{I}& K_{\ast }\\ K_{\ast }^{\mathrm{T}}& K_{\ast \ast }\end{array}\right])$$
其中，${\mathrm{X}}_{\ast }$表示预测输入，而$f_{\ast }$表示预测输出，$m_{\ast }=u_0({\mathrm{X}}_{\ast })$， K ∗ T = { k ( x 1 , X ∗ ) , … , k ( x n , X ∗ ) } K_*^\mathrm{T}=\{k(\mathrm{x}_1,\mathrm{X_*}), \ldots, k(\mathrm{x}_n,\mathrm{X}_*)\} K∗T​={k(x1​,X∗​),…,k(xn​,X∗​)},$K_{\ast \ast }=k({\mathrm{X}}_{\ast },{\mathrm{X}}_{\ast })$
我们有：
$$p(f_{\ast }|{\mathrm{X}}_{\ast },\mathrm{X},y,{\sigma }^2)=\frac{p(f_{\ast },y|{\mathrm{X}}_{\ast },\mathrm{X},{\sigma }^2)}{p(y|\mathrm{X},{\sigma }^2)}$$
所以，同样的，我们只需要考虑分子关于$f_* $的核就行了。
p ( f ∗ , y ∣ X ∗ , X , σ 2 ) ∝ e x p { − ( f ∗ − m ∗ ) T C ( f ∗ − m ∗ ) + 2 ( f ∗ − m ∗ ) T B T ( y − m ) 2 } p(f_*,y|\mathrm{X}_*,\mathrm{X}, \sigma^2) \propto \mathrm{exp}\{-\frac{(f_*-m_*)^{\mathrm{T}}C(f_*-m_*)+2(f_*-m_*)^{\mathrm{T}}B^{\mathrm{T}}(y-m)}{2}\} p(f∗​,y∣X∗​,X,σ2)∝exp{−2(f∗​−m∗​)TC(f∗​−m∗​)+2(f∗​−m∗​)TBT(y−m)​}
其中:
$$\left[\begin{array}{cc}A& B\\ B^{\mathrm{T}}& C\end{array}\right]={\left[\begin{array}{cc}K+{\sigma }^2\mathrm{I}& K_{\ast }\\ K_{\ast }^{\mathrm{T}}& K_{\ast \ast }\end{array}\right]}^{-1}$$
容易推得$f_{\ast }$的均值${\mu }_{\ast }$和协方差${\sigma }_{\ast }^2$为：
$$\begin{array}{l}{\mu }_{\ast }=-C^{-1}{B}^{\mathrm{T}}(y-m)+m_{\ast }\\ {\sigma }_{\ast }^2=C^{-1}\end{array}$$

根据Schur补和分块矩阵求逆（凸优化-p622）的性质，我们可以得到：
$$\begin{array}{c}{C}^{-1}=K_{\ast \ast }-K_{\ast }^{\mathrm{T}}(K+{\sigma }^2\mathrm{I}{)}^{-1}{K}_{\ast }\\ B^{\mathrm{T}}=-C{K}_{\ast }^{\mathrm{T}}(K+{\sigma }^2\mathrm{I}{)}^{-1}\end{array}$$
所以，我们得到：
$$\begin{array}{c}{\mu }_{\ast }({\mathrm{X}}_{\ast })=m_{\ast }+K_{\ast }^{\mathrm{T}}(K+{\sigma }^2\mathrm{I}{)}^{-1}(y-m)\\ {\sigma }_{\ast }^2=K_{\ast \ast }-K_{\ast }^{\mathrm{T}}(K+{\sigma }^2\mathrm{I}{)}^{-1}{K}_{\ast }\end{array}$$
证毕.