function hkrr_demo
%
% HKRR_DEMO.M - Heteroscedastic kernel ridge regression demo
%

%
% File        : hkrr_demo.m
%
% Date        : Monday 25h July 2005
%
% Author      : Dr Gavin C. Cawley
%
% Description : Self-contained demonstration of the heteroscedastic kernel
%               ridge regression algorithm, using Silverman's motorcycle
%               benchmark dataset.
%
% History     : 25/07/2005 - v1.00
%
% Copyright   : (c) Dr Gavin C. Cawley, July 2005.
%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
%

% This is Silverman's motorcycle benchmark dataset

data = [  2.4    0.0 ;  2.6   -1.3 ;  3.2   -2.7 ;  3.6    0.0 ;
          4.0   -2.7 ;  6.2   -2.7 ;  6.6   -2.7 ;  6.8   -1.3 ;
          7.8   -2.7 ;  8.2   -2.7 ;  8.8   -1.3 ;  8.8   -2.7 ;
          9.6   -2.7 ; 10.0   -2.7 ; 10.2   -5.4 ; 10.6   -2.7 ;
         11.0   -5.4 ; 11.4    0.0 ; 13.2   -2.7 ; 13.6   -2.7 ;
         13.8    0.0 ; 14.6  -13.3 ; 14.6   -5.4 ; 14.6   -5.4 ;
         14.6   -9.3 ; 14.6  -16.0 ; 14.6  -22.8 ; 14.8   -2.7 ;
         15.4  -22.8 ; 15.4  -32.1 ; 15.4  -53.5 ; 15.4  -54.9 ;
         15.6  -40.2 ; 15.6  -21.5 ; 15.8  -21.5 ; 15.8  -50.8 ;
         16.0  -42.9 ; 16.0  -26.8 ; 16.2  -21.5 ; 16.2  -50.8 ;
         16.2  -61.7 ; 16.4   -5.4 ; 16.4  -80.4 ; 16.6  -59.0 ;
         16.8  -71.0 ; 16.8  -91.1 ; 16.8  -77.7 ; 17.6  -37.5 ;
         17.6  -85.6 ; 17.6 -123.1 ; 17.6 -101.9 ; 17.8  -99.1 ;
         17.8 -104.4 ; 18.6 -112.5 ; 18.6  -50.8 ; 19.2 -123.1 ;
         19.4  -85.6 ; 19.4  -72.3 ; 19.6 -127.2 ; 20.2 -123.1 ;
         20.4 -117.9 ; 21.2 -134.0 ; 21.4 -101.9 ; 21.8 -108.4 ;
         22.0 -123.1 ; 23.2 -123.1 ; 23.4 -128.5 ; 24.0 -112.5 ;
         24.2  -95.1 ; 24.2  -81.8 ; 24.6  -53.5 ; 25.0  -64.4 ;
         25.0  -57.6 ; 25.4  -72.3 ; 25.4  -44.3 ; 25.6  -26.8 ;
         26.0   -5.4 ; 26.2 -107.1 ; 26.2  -21.5 ; 26.4  -65.6 ;
         27.0  -16.0 ; 27.2  -45.6 ; 27.2  -24.2 ; 27.2    9.5 ;
         27.6    4.0 ; 28.2   12.0 ; 28.4  -21.5 ; 28.4   37.5 ;
         28.6   46.9 ; 29.4  -17.4 ; 30.2   36.2 ; 31.0   75.0 ;
         31.2    8.1 ; 32.0   54.9 ; 32.0   48.2 ; 32.8   46.9 ;
         33.4   16.0 ; 33.8   45.6 ; 34.4    1.3 ; 34.8   75.0 ;
         35.2  -16.0 ; 35.2  -54.9 ; 35.4   69.6 ; 35.6   34.8 ;
         35.6   32.1 ; 36.2  -37.5 ; 36.2   22.8 ; 38.0   46.9 ;
         38.0   10.7 ; 39.2    5.4 ; 39.4   -1.3 ; 40.0  -21.5 ;
         40.4  -13.3 ; 41.6   30.8 ; 41.6  -10.7 ; 42.4   29.4 ;
         42.8    0.0 ; 42.8  -10.7 ; 43.0   14.7 ; 44.0   -1.3 ;
         44.4    0.0 ; 45.0   10.7 ; 46.6   10.7 ; 47.8  -26.8 ;
         47.8  -14.7 ; 48.8  -13.3 ; 50.6    0.0 ; 52.0   10.7 ;
         53.2  -14.7 ; 55.0   -2.7 ; 55.0   10.7 ; 55.4   -2.7 ;
         57.6   10.7 ]; 

x      = data(:,1);
t      = data(:,2);
scale  = std(t);
t      = t/scale;
ntp    = size(x,1);
X      = (min(x):0.1:max(x))';

% kernel ridge regression model

lambda = 6;
gamma  = 1;

[alpha,b,sigma] = krr(lambda, gamma, x, t);

mu    = scale*(rbf(lambda, X, x)*alpha + b);
sigma = scale*sigma;

figure(1);
clf;
plot(x, scale*t, 'go', X, mu, 'r-', X, mu+sigma, 'b--', X, mu-sigma, 'b--');
title('Kernel Rdige Regression (KRR)');
xlabel('Time (ms)');
ylabel('Acceleration (g)');

% heteroscedastic kernel ridge regression model

lambda_mu    = 8;
gamma_mu     = 1e-1;
lambda_sigma = 7;
gamma_sigma  = 1;

[alpha_mu,b_mu,alpha_sigma,b_sigma] = hkrr(lambda_mu,    ...
                                           gamma_mu,     ...
                                           lambda_sigma, ...
                                           gamma_sigma,  ...
                                           x,            ...
                                           t);

mu    = scale*(rbf(lambda_mu, X, x)*alpha_mu + b_mu);
sigma = scale*exp(rbf(lambda_sigma, X, x)*alpha_sigma + b_sigma);

figure(2);
clf;
plot(x, scale*t, 'go', X, mu, 'r-', X, mu+sigma, 'b--', X, mu-sigma, 'b--');
title('Heteroscedastic Kernel Rdige Regression (HKRR)');
xlabel('Time (ms)');
ylabel('Acceleration (g)');

% leave-one-out heteroscedastic kernel ridge regression model

lambda_mu    = 8;
gamma_mu     = 1e-1;
lambda_sigma = 7;
gamma_sigma  = 1;

[alpha_mu,b_mu,alpha_sigma,b_sigma] = loohkrr(lambda_mu,    ...
                                              gamma_mu,     ...
                                              lambda_sigma, ...
                                              gamma_sigma,  ...
                                              x,            ...
                                              t);

mu    = scale*(rbf(lambda_mu, X, x)*alpha_mu + b_mu);
sigma = scale*exp(rbf(lambda_sigma, X, x)*alpha_sigma + b_sigma);

figure(3);
clf;
plot(x, scale*t, 'go', X, mu, 'r-', X, mu+sigma, 'b--', X, mu-sigma, 'b--');
title('Leave-one-out Heteroscedastic Kernel Rdige Regression (LOOHKRR)');
xlabel('Time (ms)');
ylabel('Acceleration (g)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
% SUBFUNCTION : KRR - homoscedastic kernel ridge regression                   %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [alpha,b,sigma] = krr(lambda, gamma, x, t)

K     = rbf(lambda, x, x);
one   = ones(size(t));
R     = chol(K + gamma*eye(length(t)));
omega = R\(R'\one);
nu    = R\(R'\t);
delta = 1/max(sum(omega),realmin);
b     = sum(nu)*delta;
alpha = nu - omega*b;
sigma = std(t - K*alpha - b);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
% SUBFUNCTION : HKRR - heteroscedastic kernel ridge regression                %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [alpha_mu,b_mu,alpha_sigma,b_sigma] = hkrr(lambda_mu,gamma_mu,lambda_sigma,gamma_sigma, x, t)

% initialisation

K_mu    = rbf(lambda_mu, x, x);
K_sigma = rbf(lambda_sigma, x, x);
z_sigma = log(0.1)*ones(size(t));
sigma   = 0.1*ones(size(t));
one     = ones(size(t));

for i=1:20

   % fit model of conditional mean

   zeta     = 2*sigma.^2;
   R        = chol(K_mu + diag(gamma_mu.*zeta));
   omega    = R\(R'\one);
   nu       = R\(R'\t);
   delta    = 1/max(sum(omega),realmin);
   b_mu     = sum(nu)*delta;
   alpha_mu = nu - omega*b_mu;

   % estimate conditional mean

   mu = K_mu*alpha_mu + b_mu;

   % fit model of conditional standard deviation

   xi          = 0.5*(mu - t).^2 + realmin; 
   beta        = 4*xi.*exp(-2*z_sigma);
   eta         = z_sigma - 1./beta + 0.5;
   R           = chol(K_sigma + diag(gamma_sigma./beta));
   omega       = R\(R'\one);
   nu          = R\(R'\eta);
   delta       = 1/max(sum(omega),realmin);
   b_sigma     = sum(nu)*delta;
   alpha_sigma = nu - omega*b_sigma;
   
   % estimate conditional standard deviation

   z_sigma = K_sigma*alpha_sigma + b_sigma;
   sigma   = exp(z_sigma);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
% SUBFUNCTION : LOOHKRR - l-o-o heteroscedastic kernel ridge regression       %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [alpha_mu,b_mu,alpha_sigma,b_sigma] = loohkrr(lambda_mu,gamma_mu,lambda_sigma,gamma_sigma, x, t)

% initialisation

K_mu    = rbf(lambda_mu, x, x);
K_sigma = rbf(lambda_sigma, x, x);
z_sigma = log(0.1)*ones(size(t));
sigma   = 0.1*ones(size(t));
one     = ones(size(t));

for i=1:20

   % fit model of conditional mean

   zeta     = 2*sigma.^2;
   R        = chol(K_mu + diag(gamma_mu.*zeta));
   omega    = R\(R'\one);
   nu       = R\(R'\t);
   delta    = 1/max(sum(omega),realmin);
   b_mu     = sum(nu)*delta;
   alpha_mu = nu - omega*b_mu;

   % compute leave-one-out estimate of conditional mean

   Ri = inv(R);
   ci = sum(Ri.^2,2) - delta*omega.^2;
   mu = t - alpha_mu./max(ci,realmin);

   % fit model of conditional standard deviation

   xi          = 0.5*(mu - t).^2 + realmin; 
   beta        = 4*xi.*exp(-2*z_sigma);
   eta         = z_sigma - 1./beta + 0.5;
   R           = chol(K_sigma + diag(gamma_sigma./beta));
   omega       = R\(R'\one);
   nu          = R\(R'\eta);
   delta       = 1/max(sum(omega),realmin);
   b_sigma     = sum(nu)*delta;
   alpha_sigma = nu - omega*b_sigma;
   
   % compute conditional standard deviation

   z_sigma = K_sigma*alpha_sigma + b_sigma;
   sigma   = exp(z_sigma);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
% SUBFUNCTION : RBF - evaluate radial basis kernel function                   %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function K = rbf(lambda, x1, x2)
%
% RBF - evaluate Gaussian radial basis function kernel
%

D = sum(x1.^2,2)*ones(1,size(x2,1))+ones(size(x1,1),1)*sum(x2.^2,2)'-2*x1*x2';
K = exp(-D/(lambda^2));

% bye bye...

