高斯过程的matlab程序实现及其参数优化

高斯过程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，求最大似然估计的相反数，优化时求minimize
    outputArg1=0.5*inputy'*alpha+sum(log(diag(L)))+0.5*m*(2*pi);
    %gradient
    outputArg2 =  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==4
    n=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)';
    %方差减去sn2
    V = V-exp(2*inputloghyper(3));
    outputArg2=V;
end
end

结果：

nlZ =
   48.6337
dnlZ =

  -14.4962
    9.4411
    0.5962
mu =
   -0.9273
   -0.6795
   -0.2167
    0.2895
    0.7027
    0.9465
    0.9701
    0.7551
    0.3491
   -0.1420
   -0.5949
   -0.9024
   -0.9908
   -0.8364
   -0.4765
   -0.0000
    0.4765
    0.8364
    0.9908
    0.9024
    0.5949
    0.1420
   -0.3491
   -0.7551
   -0.9701
   -0.9465
   -0.7027
   -0.2895
    0.2167
    0.6795
    0.9273
S2 =
    0.0098
    0.0218
    0.0097
    0.0155
    0.0096
    0.0143
    0.0096
    0.0140
    0.0095
    0.0139
    0.0095
    0.0139
    0.0095
    0.0139
    0.0095
    0.0139
    0.0095
    0.0139
    0.0095
    0.0139
    0.0095
    0.0139
    0.0095
    0.0140
    0.0096
    0.0143
    0.0096
    0.0155
    0.0097
    0.0218
    0.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)

结果