机器学习&数据挖掘笔记_10（高斯过程简单理解）

　　前言：

　　GP（高斯过程）是一种自然界中普遍存在且重要的随机过程，也叫正态随机过程，在ML等领域应用比较广泛。本次实验目的是简单理解下GP，特别是要体验到GP的一个sample不再是一个普通的点，而是一个函数。实验部分完成了常见的GP的一维和二维sample的显示，常见的GP有线性GP，布朗运动，指数GP，Ornstein-Uhlenbeck过程，对称和周期性的GP。参考资料全为视频http://www.youtube.com/playlist?list=PLD0Z06AA0D2E8ZZBA中相关的部分。

　　实验基础：

　　首先来看看GP的定义：

 　　

　　GP是高斯函数的概率分布，给定S集合中的任意个自变量Si, Sj ,… Sk, 满足z(Si), z(Sj), …,z(Sk)是高斯随机变量，且z(Si), z(Sj), …,z(Sk)是多维高斯函数（可以有任意多个），记为Z，此时称随机变量集合z(Si), z(Sj), …,z(Sk)为S集合上的高斯过程。

　　如果S集合中的元素个数是有限的，则Z是否为GP可以通过穷举判断其是否为多维高斯函数。如果S集合中元素的个数是无限的，则Z不能通过穷举获得，但是如果Z中变量和某个高斯变量有直接联系的话，则Z也有可能是S上的一个GP。

　　下面来看看GP的存在性定理以及该怎样构造 GP：

 　　

　　存在定理说明，对任意集合S中的单个元素都存在某个均值函数，以及对任意集合S中的2个元素都存在某个核函数（即协方差函数），则在S上一定存在一个高斯过程Z(t)，其元素具有类似S形式的均值和方差。所以在给定集合S后，我们只需要给出一个一元的均值函数，一个二元的核函数表达式，就可以构造出一个高斯过程。常见的高斯过程有：

　　Random planes, Brownian motion, squared exponential GP, Ornstein-Uhlenbeck, a periodic GP, and a symmetric GP.

　　它们的均值函数和协方差函数分别定义如下：

 　　　　

　　

　　这里讲一下通常情况对高斯函数采样的方法：首先我们知道任何一个高斯函数都可以写成标准高斯函数的线性组合，因此只要能够对标准高斯函数进行采样就OK了。其方法为：计算出了标准高斯函数的分布函数，用[0,1]均匀分布随机发生器选择随机的值y，当做标准高斯函数的函数值，然后找到分布函数下对应的S就可以了，该点即为我们所需要的sample。

　　程序中通过svd分解后的变量来重构z的理论暂时没有理解。

　　实验结果：

　　首先看看1d情况下的GP sample，其线性GP如下（图像表明，符合该条件的GP中的一个sample是一条直线）：

 　　

　　布朗运动的结果如下（有点随机游动的感觉）：

 　　

　　指数GP结果为（特点是很平滑）：

 　　

　　Ornstein-Uhlenbeck过程如下（和布朗运动类似）：

 　　

　　周期性的GP：

 　　

　　对称性的GP：

 　　

　　下面是2d情况下的结果。线性的GP(是一个平面)：

 　　

　　squared exponential：

 　　

　　Ornstein-Uhlenbeck过程：

 　　

　　2d的处理耗时和耗内存比较大，所以这里将显示的样本点个数减小了。

　　实验主要部分代码及注释：

GP_1d.m：

% GP_1d.m

% 选择核函数(协方差函数)，其均值函数默认为常量0.
kernel = 6;
switch kernel
    case 1; k = @(x,y) 1*x'*y; %Linear
    case 2; k = @(x,y) 1*min(x,y); % Brownian  
    case 3; k = @(x,y) exp(-100*(x-y)'*(x-y)); %squared exponential
    case 4; k = @(x,y) exp(-1*sqrt((x-y)'*(x-y))); %Ornstein-Uhlenbeck
    case 5; k = @(x,y) exp(-1*sin(5*pi*(x-y)).^2); %A periodic GP
    case 6; k = @(x,y) exp(-100*min(abs(x-y),abs(x+y)).^2); %A symmetric GP
end

% 选择需显示的点x，即集合S的一部分
x = (-1:0.005:1);
n = length(x);

% 构造协方差矩阵
C = zeros(n,n);
for i = 1:n
    for j = 1:n
        C(i,j) = k(x(i),x(j));
    end
end

% 对GP进行采样
rn = randn(n,1);%产生n个0~1之间的随机数,满足正态分布
[u,s,v] = svd(C); %svd分解rn矩阵，s为奇异值矩阵，u为奇异向量.C=usv'
z = u*sqrt(s)*rn; %z为什么这么表示,理论是？？

% 画出GP的一个sample
figure(1);hold on; clf
plot(x,z,'.-');
% axis([0,1,-2,2]);

GP_2d.m：

% GP_2d.m

% 选择核函数(协方差函数)，其均值函数默认为常量0.
kernel = 1;
switch kernel
    case 1; k = @(x,y) 1*x'*y; %Linear
    case 2; k = @(x,y) exp(-100*(x-y)'*(x-y)); %squared exponential
    case 3; k = @(x,y) exp(-1*sqrt((x-y)'*(x-y))); %Ornstein-Uhlenbeck
end

% 选择需显示的点points，二维的
points = (0:0.02:1)';
[U,V] = meshgrid(points,points);
x = [U(:) V(:)]';
n = size(x,2);

% 构造协方差矩阵
C = zeros(n,n);
for i = 1:n
    for j = 1:n
        C(i,j) = k(x(:,i),x(:,j));
    end
end

% 对GP进行采样
rn = randn(n,1);%产生n个0~1之间的随机数,满足正态分布
[u,s,v] = svd(C); %svd分解rn矩阵，s为奇异值矩阵，u为奇异向量.C=usv'
z = u*sqrt(s)*rn; %z为什么这么表示,理论是？？

% 画出GP的一个sample
figure(2); clf
Z = reshape(z,sqrt(n),sqrt(n));
surf(U,V,Z);

　　参考资料：

     http://www.youtube.com/playlist?list=PLD0Z06AA0D2E8ZZBA