高斯过程matlab编程实战详解

本次结论公式均来自《Gaussian Processes for Machine Learning》-> http://www.gaussianprocess.org/gpml/

P19

P114

OK,开始吧。

问题记录：

1、如果某些情况下出现核矩阵不对称，可通过(H+H')/2进行处理

2、矩阵求逆不正定，加一个正则化项，也叫吉洪诺夫正则化，就是单位矩阵eye(n)*epsilon，但是epsilon需要特别小，太大就变成回归形式了，但是多小是个问题，1e-7/（1e-22*n）都有，以后有时间一定要好啊后研究研究，感觉写程序这个东西很难受

样本：

x=(-7.5:1:7.5)';
y=sin(x);

预测点：

xstar=(-7.5:0.5:7.5)';

超参数初始值：

loghyper = [log(1.0); log(1.0);  log(0.1)];  %分别为length-scale，sigmaf，sigman并取ln对数

求最大似然估计和倒数或者预测均值和方差

[nlZ dnlZ] =  GPR(loghyper, x,y)
[mu, S2] = GPR(loghyper,  x,y,xstar)

协方差中的（欧式距离）^2函数

function C =  SQDIST(Sample,Pre);
%求样本的欧氏距离并放入对应的协方差矩阵中，可以看作
% for i=1:n
%     for j=1:n
%         K(i,j)=norm(x(i,:)-x(j,:))^2;
%     end
% end
if nargin==1[n,D1] = size(Sample);C=zeros(n);for d = 1:D1    C=C+(repmat(Sample(:,d),  1 ,m) - repmat(Sample(:,d)',  n,1)).^2;   end
end
%求样本与样本点间的欧氏距离并放入对应的协方差矩阵中,可以看作
% for i=1:m
%     for j=1:n
%         K(i,j)=norm(x(i,:)-xpre(j,:))^2;
%     end
% end
if nargin==2[n,D1] = size(Sample);[m,D2] = size(Pre);C=zeros(n,m);for d = 1:D1    C=C+(repmat(Sample(:,d),  1 ,m) - repmat(Pre(:,d)',  n,1)).^2;   end
end
end

函数->求预测值和方差或者求最大似然估计及其梯度

function  [outputArg1,outputArg2] =  GPR(inputloghyper,  inputx,inputy,inputxstar)
% m个dim维的样本
% inputx:m*dim
% inputy:m*1
% inputxstar:n*dim
%  inputloghyper:log(hyper)=[theta(1),theta(2),theta(3)]
ell = exp(inputloghyper(1));                             % characteristic length-scale
sf2 = exp(2*inputloghyper(2));                           % sigmaf
sn2 = exp(2*inputloghyper(3));                           % sigman:noise variance
[m,dim]=size(inputx);
%按照公式2.31求协方差矩阵%
Kxx=zeros(m,m);
KKxx=sf2*exp(-0.5*SQDIST(inputx/ell,inputx/ell));
Kxx=Kxx+KKxx;
KKxx=sn2*eye(size(inputx,1));
Kxx=Kxx+KKxx;
%cholesky分解
L=chol(Kxx)';
alpha=L'\(L\inputy);
%求最大似然估计及其梯度
if nargin==3% minlogpyx，求最大似然估计的相反数，优化时求minimizeoutputArg1=0.5*inputy'*alpha+sum(log(diag(L)))+0.5*m*(2*pi);%gradientoutputArg2 =  zeros(size(inputloghyper));W =  L'\(L\eye(m))-alpha*alpha';%求K的关于超参数的偏导数，这里注意，我们的参数是lnx，这里sf2的偏导数是exp(2*logtheta)' =2*exp(2*logtheta)=2*sf2，而不是sf^2'=2*sf=2*exp(logtheta)dell =  sf2*exp(-SQDIST(inputx/ell)/2).*SQDIST(inputx/ell);  dsf =  2*sf2*exp(-SQDIST(inputx/ell)/2);dsn =  2*sn2*eye(size(inputx,1));%求迹，这里参考GPML，但是为什么这样算还不知道，直接求迹不对，望大神指点迷津outputArg2(1)=sum(sum(W.*dell))/2;outputArg2(2)=sum(sum(W.*dsf))/2;outputArg2(3)=sum(sum(W.*dsn))/2;
end
%求预测值和方差
if nargin==4n=size(inputxstar,1);Kxxstar=zeros(m,n);KKxxstar=sf2*exp(-0.5*SQDIST(inputx/ell,inputxstar/ell));Kxxstar=Kxxstar+KKxxstar;%求均值fmean=Kxxstar'*alpha;outputArg1=fmean;v = L\Kxxstar;V=sf2+sn2-sum(v.*v)';%方差减去sn2V = V-exp(2*inputloghyper(3));outputArg2=V;
end
end

结果：

nlZ =48.6337
dnlZ =-14.49629.44110.5962
mu =-0.9273-0.6795-0.21670.28950.70270.94650.97010.75510.3491-0.1420-0.5949-0.9024-0.9908-0.8364-0.4765-0.00000.47650.83640.99080.90240.59490.1420-0.3491-0.7551-0.9701-0.9465-0.7027-0.28950.21670.67950.9273
S2 =0.00980.02180.00970.01550.00960.01430.00960.01400.00950.01390.00950.01390.00950.01390.00950.01390.00950.01390.00950.01390.00950.01390.00950.01400.00960.01430.00960.01550.00970.02180.0098

画图：mu+-2sigma

figure
f = [mu+2*sqrt(S2);flipdim(mu-2*sqrt(S2),1)];
fill([xstar; flipdim(xstar,1)], f, [7 7 7]/8, 'EdgeColor', [7 7 7]/8);
hold on
plot(xstar,mu,'k-','LineWidth',2);
plot(x, y, 'k+', 'MarkerSize', 17);
figure
errorbar(xstar, mu, 2*sqrt(S2), 'g');
hold on
plot(x, y, 'k+', 'MarkerSize', 17)

结果