高斯过程回归

2019-09-10 15:38 作者:露保协 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)是统计文章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.中的两个理论支柱。关于生化反应系统（Markov跳过程）的Bayesian inference在Wilkinson D J. Stochastic modelling for systems biology[M]. CRC press, 2011.一书的最后两章有讲（这本书前面部分是概率模型，虽然对我而言都很熟知了，但是初学者还是可以看的），Gaussian process regression则在这片笔记中总结。GP这东西在机器学习中也很有用，以后也可能会做这方面的工作。不过现在这篇笔记纯粹是为了解决细胞生物问题，算是一股清流。

概述

Gaussian process regression是一个exploit unknown functions的手段。从原理上来说，它是一个非参数的Bayes推断，在函数空间上（未知函数视为随机元），由先验分布加入信息得到后验分布。函数空间上的分布即随机过程，因为是无限维的，所以是非参数的。这样一个分布（测度）在本文中均考虑为Gauss过程（或Gauss随机场），它有很多好性质，最主要是好算。Gauss过程其实就是个无穷维的多元正态分布，可以用期望函数和协方差函数（kernel）完全刻画；只不过此时要由Kolmogorov's extension theorem保证了。

我们有一个complex model（可以是一个computer model，比如是一个完整的Gillespie算法），它是一个黑箱，输入一个input，给出一个output。我们只知道有这样一个函数关系f(x)，以及它在特定点上的值，但是因为不知道黑箱内部是什么，所以无法完整地知道f(.)。对于这个未知函数，我们有一个先验的猜测（prior distribution，取决于问题，比如输入相差小，输出也相差小，那么kernel应该是随|t-t'|衰减的），在加入信息修正之后，变成一个更精确的posterior distribution。

我算是贝叶斯学派的，碰见未知的东西别的不说先假设一个先验分布，下文也全部按照这一思想来。

simulator与emulator的区别

simulator就是前面说的computer model，它用精确的黑箱机制输出output。但是我们在想：

不知道f(.)的具体形式，只能数值计算；
更重要的是，simulator可能需要很长很长时间来算，比如一个长时间的Gillespie。但是在Bayesian inference中，我们用Metropolis-Hastings algorithm采样posterior ditribution时，每一步都需要一个Gillespie，复杂度O(MN)，无法接受。

所以：我们想要一个simulator的statistical approximation，它能够：

从simulator的统计数据中训练出来；
给出f(.)的一个Bayesian的猜测，这个猜测本身也是统计的；也就是说，它包含了f(.)的posterior distibution。所以我们可以知道误差多大，95%置信区间在哪里，等等[Lagrangian interpolation无法做到的]事情。

这个statistical approximation就叫做emulator。

中文翻译里规定simulate叫模拟，emulate叫仿真，很贴切。

例子：果蝇胚胎发育中的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.），本质上讲的就是一句话：果蝇胚胎上的细胞通过蛋白浓度找自己的位置，用的是Bayesian方法（虽然作者似乎没有弄清这个重点）。

对于这个例子，simulator就是黑箱内部具体的分子机制，emulator就是人为构建的Bayesian maximun maximum a posteriori算法，emulator能够在不知道分子机制的情况下近似仿真simulator。而人为构建的Bayesian的emulator能够好地近似simulator，就说明这个simulator本身也是近似Bayesian的，这大概是进化压力的作用。生命很多地方体现出的Bayesian都很神奇，比如狼来了这种例子，再比如最近关于猕猴身体感知的文章（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.，有空读一下），应该都是进化出来有利于生存的。

GPR的本质性在于哪儿？

平时说的线形回归，或者各种插值，某种特殊形式函数的回归，都是「有限维」的，即有有限的几个参数构成的f(.)备选函数类。这当然是基于对黑箱的认知，如果它的机制是线性+噪声，拿当然很好。但是如果没有备选函数类的认知，那么要在全函数空间上找到一个函数，这就是一个「无限维」的问题。GPR本质就在于能解决无限维（non-parametric）的问题。

比如说，我们想知道一个细胞中某种反馈的函数形式，一般也就参数化地假设成Hill函数。但是如果没有Hill函数的道理的话，我们总不可能用拉格朗日插值。我们要探索的就是「函数形式」本身，此时只能用非参数的统计。

在GPR中，我们假设输入和输出的关系为

（分布用mathcal写显得有逼格）

统计学中的一个小trick

我们不太喜欢supported on a certain subset of R的函数或者分布。为此可以做一些变换，比如指数变换（对数正态分布），logit变换（(0,1)->R）。在supported on R的函数上面就可以随便做了，正态也好，不用关心它会不会掉出去也罢。

不同kernel对Gaussian process的影响

不同kernel表示不同时间点上GP两个点的关联程度。如果它随|t-t'|衰减，就有可能产生光滑曲线；如果是min(t,t')，那就是完全不光滑的Brownian motion。实际上对于一个computer model，应该满足输入相差越小，输出关联性/informative connection越大这个基本假设，这就对kernel的选取提出了要求。

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

首先，问题在于怎么从一个函数空间上采样呢？我们实际上永远无法在计算机上做无穷维的事情；我们只能离散地取很多点，计算这些点上的函数值——用多维Gauss分布的采样（先从GP计算出均值向量与协方差矩阵）。比如可以用以下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有什么区别，以及应该选取哪个作为prior。

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

也就是我们熟悉的样子。

一个stationary kernel定义为x-x'的函数，它产生stationary的GP。BM当然不是stationary的。近一步如果是|x-x'|的函数则称为isotropic，此时也叫radial basis functions (RBFs)。

上面说的是一般的函数理论里的kernel。但是对于covariace function，当然要求kernel具有对称性，即k(x,x')=k(x',x)，所以必须是isotropic的。

如果是x\cdot x'的函数，则称为dot product kernel，它是旋转不变的，也是对称的。

给定一串输入x_i，K_{ij}=k(x_i,x_j)称为Gram matrix，它是一个协方差矩阵，当然这要求kernel是半正定的（Mercer定理）。

下面只考虑stationary kernel。记\tau=x-x'，r=|\tau|。对于k(\tau)，它有一个（非负）的（功率）谱密度S(s)（Bochner’s theorem）:

谱密度提供了函数光滑性的信息。

squared exponential (SE) kernel的形式为：

其中的l和sigma_f称为hyperparameter。sigma_f无所谓的，之后就略去了。l称为特征长度，即表示关联的长度尺度。它的功率谱也是同样形式的平方指数衰减。Stein认为这个kernel太光滑了，不适合很多物理体系，Matern kernel才更好；然而大家都这么用（SE kernel是用的最多的）。

SE kernel的采样如下：