貝葉斯優化

目錄

• 簡介
  • 貝葉斯優化框架
• 概率代理模型
  • 參數模型
    • 湯普森采樣和Beta-Bernouli模型
  • 線性模型(Linear models)
• 非參數模型
  • 高斯過程
  • 常用的一些kernels
  • 邊際似然
  • 復雜度
  • sparse pseudo-input Gaussian processes (SPGP)
  • Sparse spectrum Gaussian processes(SSGP)
  • 隨機森林
• 采集函數
  • 基於提升的策略
    • PI (probability of improvement)
    • EI(expected improvement)
  • UCB(upper confidence bound)
  • 基於信息的策略
  • 采集函數的組合
• 出於實際的考量
  • 處理超參數
  • 采集函數的優化
• 數值實驗
  • Beta-Bernoulli
  • 高斯過程
• 代碼
• 附錄
  • Beta-Bernoulli 模型的推導
  • 線性模型 后驗分布$p(\mathrm{w},\sigma^2|D_n)$的推導
  • $y$邊際分布推導
  • $f_*$的后驗分布

[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^{*}}=\mathop{argmax}\limits_{\mathrm{x}\in \mathcal{X}}f(\mathrm{x}) \]

注: 使用"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(\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 \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}\)的先驗分布，因為這是其共軛先驗分布。

\[p(\mathrm{w}|\alpha, \beta)=\prod \limits_{a=1}^{K}\mathrm{Beta}(w_a|\alpha, \beta) \]

定義：

\[\begin{array}{ll} n_{a,0}=\sum \limits_{i=1}^{n}\mathbb{I}(y_i=0, a_i=a)\\ n_{a,1}=\sum \limits_{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 \limits_{a=1}^{K} \mathrm{Beta}(w_a|\alpha+n_{a,1}, \beta+n_{a,0}) \]

上述推導見附錄。
從上述也能發現，超參數\(\alpha, \beta\)表示的治愈數和失敗數。下圖是以\(\mathrm{Beta}(w|2,2)\)為先驗的一個例子。

湯普森采樣-wiki
那么在\(D_n\)的基礎上，如果找\(a_{n+1}\)呢。以下采用的是湯普森采樣（或后驗采樣）：

\[\alpha_{n+1}=\mathop{argmax}\limits_{a} f_{\widetilde{\mathrm{w}}}(a) \]

\(\widetilde{\mathrm{w}}\sim p(\mathrm{w}|D_n)\),即\(\widetilde{\mathrm{w}}\)從\(\mathrm{w}\)的后驗分布中采樣得到。
該模型的好處是：

• 只有\(\alpha,\beta\)倆個超參數
• 是對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}\)成了權重向量。這只是代理模型的一部分，因為並沒有體現“概率”的部分。

\[y_a \sim N(\mathrm{x}_a^{\mathrm{T}}\mathrm{w}, \sigma^2) \]

組合\(a\)的觀測值為\(y_a\)，服從正態分布。很自然地，我們同樣選擇共軛先驗分布作為\(\mathrm{w}, \sigma^2\)的先驗分布：normal_inverse_gamma-wiki

\[\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} \]

\(\mathrm{NIG}\)分布有4個超參數，而\(\mathrm{w}, \sigma^2\)的后驗分布（\(D_n\)的條件下）滿足：

其中\(\mathrm{X}\)的第\(i\)行為\(\mathrm{x}_{\alpha_i}\)。
推導見附錄。

關於\(\alpha_{n+1}\)的選擇，同樣可以采用湯普森采樣：

\[\alpha_{n+1}=\mathop{argmax}\limits_{a} \mathrm{x}_a^{\mathrm{T}}\widetilde{\mathrm{w}} \]

其中\(\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})\)常常為：

\[\mathrm{exp}\{ -\frac{(\mathrm{x-z_j})^{\mathrm{T}}\Lambda(\mathrm{x-z_j})}{2} \} \]

和

\[\mathrm{exp}\{ -i\mathrm{x^T}w_j\} \]

這里，\(\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\)的多維正態分布。那么，\(y\sim \mathcal{N}(\mathrm{Xw},\sigma^2\mathrm{I})\)，我們可以積分掉\(\mathrm{w}\)從而獲得\(y\)的一個邊際分布:

推導見附錄。
就像先前已經提到過的，我們可以引入\(\Phi = (\phi_i,\ldots,\phi_J)^{T}\),
將資料（設計）矩陣\(X\)映射到\(\mathbb{R}^{n\times J}\)，如此一來，相應的邊際分布也需發生變化：

注意到\(\Phi V_0 \Phi^{T}\)，事實上，我們不需要特別指明\(\Phi\)，而只需通過kernel.

\(K\)將是\(\Phi V_0 \Phi^{T}\)的一個替代。采用這個策略，比原先在計算上和可解釋性上更有優勢。
不過，還有另外一個問題，如果去尋找下一個評估點呢。尋找下一個評估點，需要我們做預測，但是上面的邊際分布顯然是無法進行預測的。不過，我們可以通過條件公式來獲得：

\(\mathrm{X}_{*}\)是新的位置,而\(y_{*}\)是相應的預測，二者都可以是向量。
分子部分是一個聯合的高斯分布。到此，我們實際上完成了一個簡易的高斯過程，下面正式介紹高斯過程。

高斯過程

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

高斯過程\(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_*\)的分布。
首先，

\[p(y|\mathrm{X}, \sigma^2) = \mathcal{N}(y|m, K+\sigma^2\mathrm{I}) \]

我們並沒有給出這個證明，因為這個不難驗證。接下來，為了預測，我們需要求后驗分布。論文此時並沒有選擇\(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\)是一個對角矩陣，其對角線元素為\(\theta_i^2,i=1,\ldots,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(nm^2+m^3)\)。缺點是，參數相對比較多，容易產生“過擬合”現象。

Sparse spectrum Gaussian processes(SSGP)

由Bochner定理得，任何穩定得kernel都有一個正定得傅里葉譜表示：

之后，通過蒙特卡洛方法，采樣m個樣本頻率，來近似估計上訴的積分。從而獲得近似的協方差函數，當數據集較小時，SSGP同樣易產生“過擬合”現象。

隨機森林

隨機森林 - Poll的筆記

隨機森林可以作為高斯過程的一種替代。缺點是，數據缺少的地方，估計的並不准確（邊際更是常數）。另外，由於隨機森林不連續，也就不可微，所以無法采用梯度下降（或上升）的方法來更新參數。另外不解的是，隨機森林的參數，即便我們給了一個先驗分布，可其后驗分布如何求呢？

采集函數

首先我們有一個效用函數\(U:\mathbb{R}^d \times \mathbb{R} \times \Theta \rightarrow \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 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)

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

而相應的的采集函數是：

其中\(\phi\)是標准正態分布的概率密度函數。式子通過積分變量替換可以推得。
實際上\(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}^*\)的后驗分布\(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)d\mathrm{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{\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為 \(\theta_0^2 exp\{-\frac{1}{2}r^2 \}\)，采集函數為\(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{Beta}(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 \limits_{a=1}^{k}p(\mathrm{w}_a|D)\\ &\propto \prod \limits_{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{Beta}(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|\mathrm{w},\sigma^2)\propto \frac{1}{\sigma^n}\mathrm{exp}\{-\frac{(y-\mathrm{Xw})^{\mathrm{T}}(y-\mathrm{Xw})}{2\sigma^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} \} \]

先配\(\mathrm{w}\)，再配第二部分，即可得：

\(y\)邊際分布推導

如果我們令\(f=\mathrm{Xw}\)，那么\(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^2I) \mathcal{N} (f|0, \mathrm{XV_0X^{T}})df \end{array} \]

這個積分比先前的積分要容易求解（至少前面的那個我直接求不出來）。
依舊采用核方法：

\[\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} \]

其中:

\[\begin{array}{l} \mu = (\Sigma^{-1}+\sigma^{-2})^{-1}\sigma^{-2}y \\ \Sigma = \mathrm{XV_0X^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_*\)的后驗分布

我們不加推導地給出：

\[p(y|\mathrm{X}, \sigma^2) = \mathcal{N}(y|m, K+\sigma^2\mathrm{I}) \]

根據高斯過程的性質，存在如下的聯合分布：

\[\left [ \begin{array}{c} y\\ f_* \end{array} \right ] \sim \mathcal{N}\Big( \left [ \begin{array}{c} m\\ m_* \end{array} \right ], \left [ \begin{array}{cc} K+\sigma^2\mathrm{I} & K_*\\ K_*^\mathrm{T} & K_{**} \end{array} \right ] \Big) \]

其中，\(\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}_*)\)
我們有：

\[p(f_*|\mathrm{X}_*, \mathrm{X}, y, \sigma^2)=\frac{p(f_*,y|\mathrm{X}_*,\mathrm{X}, \sigma^2)}{p(y|\mathrm{X}, \sigma^2)} \]

所以，同樣的，我們只需要考慮分子關於$f_* $的核就行了。

\[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}\} \]

其中:

\[\left [\begin{array}{cc} A & B \\ B^{\mathrm{T}} & C \end{array} \right ]= \left [\begin{array}{cc} K+\sigma^2\mathrm{I} & K_* \\ K_*^{\mathrm{T}}& K_{**} \end{array} \right ]^{-1} \]

容易推得\(f_*\)的均值\(\mu_*\)和協方差\(\sigma_*^2\)為：

\[\begin{array}{l} \mu_* = -C^{-1}B^{\mathrm{T}}(y-m)+m_*\\ \sigma_*^2=C^{-1} \end{array} \]

根據Schur補和分塊矩陣求逆（凸優化-p622）的性質，我們可以得到：

\[\begin{array}{c} C^{-1} = K_{**}-K_*^{\mathrm{T}}(K+\sigma^2\mathrm{I})^{-1}K_*\\ B^{\mathrm{T}} = -CK_*^{\mathrm{T}}(K+\sigma^2\mathrm{I})^{-1} \end{array} \]

所以，我們得到：

\[\begin{array}{c} \mu_*(\mathrm{X}_*)=m_*+K_*^{\mathrm{T}}(K+\sigma^2\mathrm{I})^{-1}(y-m)\\ \sigma_*^2=K_{**}-K_*^{\mathrm{T}}(K+\sigma^2\mathrm{I})^{-1}K_* \end{array} \]

證畢.