高斯過程回歸

2019-09-10 15:38 作者:露保協(xié) 0人讀過 | 我要投稿

本筆記的來源

https://peterroelants.github.io/posts/gaussian-process-tutorial/的三篇文章
Schulz E, Speekenbrink M, Krause A. A tutorial on Gaussian process regression: Modelling, exploring, and exploiting functions[J]. Journal of Mathematical Psychology, 2018, 85: 1-16.
Rasmussen C E. Gaussian processes in machine learning[C]//Summer School on Machine Learning. Springer, Berlin, Heidelberg, 2003: 63-71.的第一章與第四章

Gaussian process regression（GPR）和Bayesian inference(for latent variable model)是統(tǒng)計文章Henderson D A, Boys R J, Krishnan K J, et al. Bayesian emulation and calibration of a stochastic computer model of mitochondrial DNA deletions in substantia nigra neurons[J]. Journal of the American Statistical Association, 2009, 104(485): 76-87.中的兩個理論支柱。關(guān)于生化反應(yīng)系統(tǒng)（Markov跳過程）的Bayesian inference在Wilkinson D J. Stochastic modelling for systems biology[M]. CRC press, 2011.一書的最后兩章有講（這本書前面部分是概率模型，雖然對我而言都很熟知了，但是初學(xué)者還是可以看的），Gaussian process regression則在這片筆記中總結(jié)。GP這東西在機器學(xué)習(xí)中也很有用，以后也可能會做這方面的工作。不過現(xiàn)在這篇筆記純粹是為了解決細(xì)胞生物問題，算是一股清流。

概述

Gaussian process regression是一個exploit unknown functions的手段。從原理上來說，它是一個非參數(shù)的Bayes推斷，在函數(shù)空間上（未知函數(shù)視為隨機元），由先驗分布加入信息得到后驗分布。函數(shù)空間上的分布即隨機過程，因為是無限維的，所以是非參數(shù)的。這樣一個分布（測度）在本文中均考慮為Gauss過程（或Gauss隨機場），它有很多好性質(zhì)，最主要是好算。Gauss過程其實就是個無窮維的多元正態(tài)分布，可以用期望函數(shù)和協(xié)方差函數(shù)（kernel）完全刻畫；只不過此時要由Kolmogorov's extension theorem保證了。

我們有一個complex model（可以是一個computer model，比如是一個完整的Gillespie算法），它是一個黑箱，輸入一個input，給出一個output。我們只知道有這樣一個函數(shù)關(guān)系f(x)，以及它在特定點上的值，但是因為不知道黑箱內(nèi)部是什么，所以無法完整地知道f(.)。對于這個未知函數(shù)，我們有一個先驗的猜測（prior distribution，取決于問題，比如輸入相差小，輸出也相差小，那么kernel應(yīng)該是隨|t-t'|衰減的），在加入信息修正之后，變成一個更精確的posterior distribution。

我算是貝葉斯學(xué)派的，碰見未知的東西別的不說先假設(shè)一個先驗分布，下文也全部按照這一思想來。

simulator與emulator的區(qū)別

simulator就是前面說的computer model，它用精確的黑箱機制輸出output。但是我們在想：

不知道f(.)的具體形式，只能數(shù)值計算；
更重要的是，simulator可能需要很長很長時間來算，比如一個長時間的Gillespie。但是在Bayesian inference中，我們用Metropolis-Hastings algorithm采樣posterior ditribution時，每一步都需要一個Gillespie，復(fù)雜度O(MN)，無法接受。

所以：我們想要一個simulator的statistical approximation，它能夠：

從simulator的統(tǒng)計數(shù)據(jù)中訓(xùn)練出來；
給出f(.)的一個Bayesian的猜測，這個猜測本身也是統(tǒng)計的；也就是說，它包含了f(.)的posterior distibution。所以我們可以知道誤差多大，95%置信區(qū)間在哪里，等等[Lagrangian interpolation無法做到的]事情。

這個statistical approximation就叫做emulator。

中文翻譯里規(guī)定simulate叫模擬，emulate叫仿真，很貼切。

例子：果蠅胚胎發(fā)育中的simulator與emulator

就是今年在cell上那篇文章（Petkova M D, Tka?ik G, Bialek W, et al. Optimal decoding of cellular identities in a genetic network[J]. Cell, 2019, 176(4): 844-855. e15.），本質(zhì)上講的就是一句話：果蠅胚胎上的細(xì)胞通過蛋白濃度找自己的位置，用的是Bayesian方法（雖然作者似乎沒有弄清這個重點）。

對于這個例子，simulator就是黑箱內(nèi)部具體的分子機制，emulator就是人為構(gòu)建的Bayesian maximun?maximum a posteriori算法，emulator能夠在不知道分子機制的情況下近似仿真simulator。而人為構(gòu)建的Bayesian的emulator能夠好地近似simulator，就說明這個simulator本身也是近似Bayesian的，這大概是進化壓力的作用。生命很多地方體現(xiàn)出的Bayesian都很神奇，比如狼來了這種例子，再比如最近關(guān)于獼猴身體感知的文章（Fang W, Li J, Qi G, et al. Statistical inference of body representation in the macaque brain[J]. Proceedings of the National Academy of Sciences, 2019: 201902334.，有空讀一下），應(yīng)該都是進化出來有利于生存的。

GPR的本質(zhì)性在于哪兒？

平時說的線形回歸，或者各種插值，某種特殊形式函數(shù)的回歸，都是「有限維」的，即有有限的幾個參數(shù)構(gòu)成的f(.)備選函數(shù)類。這當(dāng)然是基于對黑箱的認(rèn)知，如果它的機制是線性+噪聲，拿當(dāng)然很好。但是如果沒有備選函數(shù)類的認(rèn)知，那么要在全函數(shù)空間上找到一個函數(shù)，這就是一個「無限維」的問題。GPR本質(zhì)就在于能解決無限維（non-parametric）的問題。

比如說，我們想知道一個細(xì)胞中某種反饋的函數(shù)形式，一般也就參數(shù)化地假設(shè)成Hill函數(shù)。但是如果沒有Hill函數(shù)的道理的話，我們總不可能用拉格朗日插值。我們要探索的就是「函數(shù)形式」本身，此時只能用非參數(shù)的統(tǒng)計。

在GPR中，我們假設(shè)輸入和輸出的關(guān)系為

（分布用mathcal寫顯得有逼格）

統(tǒng)計學(xué)中的一個小trick

我們不太喜歡supported on a certain subset of R的函數(shù)或者分布。為此可以做一些變換，比如指數(shù)變換（對數(shù)正態(tài)分布），logit變換（(0,1)->R）。在supported on R的函數(shù)上面就可以隨便做了，正態(tài)也好，不用關(guān)心它會不會掉出去也罷。

不同kernel對Gaussian process的影響

不同kernel表示不同時間點上GP兩個點的關(guān)聯(lián)程度。如果它隨|t-t'|衰減，就有可能產(chǎn)生光滑曲線；如果是min(t,t')，那就是完全不光滑的Brownian motion。實際上對于一個computer model，應(yīng)該滿足輸入相差越小，輸出關(guān)聯(lián)性/informative connection越大這個基本假設(shè)，這就對kernel的選取提出了要求。

“The kernel of a Gaussian process encodes our assumptions about the function which we wish to learn.”

首先，問題在于怎么從一個函數(shù)空間上采樣呢？我們實際上永遠無法在計算機上做無窮維的事情；我們只能離散地取很多點，計算這些點上的函數(shù)值——用多維Gauss分布的采樣（先從GP計算出均值向量與協(xié)方差矩陣）。比如可以用以下matlab代碼：

%sampling from different kernels(GP on [0,1])

t=0.001:0.001:1;

mu=zeros(1,1000);

sigma=zeros(1000,1000);

for i=1:1000

? ? for j=1:i

? ? ? ? x=i/1000;

? ? ? ? y=j/1000;

? ? ? ? k=exp(-abs(x-y)^2);%set kernel

? ? ? ? sigma(i,j)=k;

? ? ? ? sigma(j,i)=k;

? ? end

end

figure(1)

hold on

z=mvnrnd(mu,sigma,5);

for i=1:5

?? plot(t,z(i,:)) %plot 5 sampled functions

end

有了這一工具，接下來就可以具體看一下不同的kernel有什么區(qū)別，以及應(yīng)該選取哪個作為prior。

首先是k(x,x')=min(x,x')，也就是布朗運動。

也就是我們熟悉的樣子。

一個stationary kernel定義為x-x'的函數(shù)，它產(chǎn)生stationary的GP。BM當(dāng)然不是stationary的。近一步如果是|x-x'|的函數(shù)則稱為isotropic，此時也叫radial basis functions?(RBFs)。

上面說的是一般的函數(shù)理論里的kernel。但是對于covariace function，當(dāng)然要求kernel具有對稱性，即k(x,x')=k(x',x)，所以必須是isotropic的。

如果是x\cdot x'的函數(shù)，則稱為dot product kernel，它是旋轉(zhuǎn)不變的，也是對稱的。

給定一串輸入x_i，K_{ij}=k(x_i,x_j)稱為Gram matrix，它是一個協(xié)方差矩陣，當(dāng)然這要求kernel是半正定的（Mercer定理）。

下面只考慮stationary kernel。記\tau=x-x'，r=|\tau|。對于k(\tau)，它有一個（非負(fù)）的（功率）譜密度S(s)（Bochner’s theorem）:

譜密度提供了函數(shù)光滑性的信息。

squared exponential?(SE) kernel的形式為：

其中的l和sigma_f稱為hyperparameter。sigma_f無所謂的，之后就略去了。l稱為特征長度，即表示關(guān)聯(lián)的長度尺度。它的功率譜也是同樣形式的平方指數(shù)衰減。Stein認(rèn)為這個kernel太光滑了，不適合很多物理體系，Matern kernel才更好；然而大家都這么用（SE kernel是用的最多的）。

SE kernel的采樣如下：