function [pr,sdev] = gp2_idencov(p,M,y,T)

% GP2_IDENCOV       Computest the prediction and the covariance matrix of
% the posterior joint probability distribution for a Gaussian process.
%
% [PR,SD] = GP2_IDENCOV(p,M,y,T)
%
% This function computes the prediction and the covariance of the posterior
% Gaussian process, given the prior information. The prior covariance 
% function used is of the form,
%
%           C(zi,zj) = a*{exp([zi-zj]'*D*[zi-zj]) + v}
%
% Normally, 95% confidence intervals corresponds to 2 times standard
% deviation, sd, for a Gaussian distribution. In this script, the 
% covariance matrix is calculated, instead of the standard deviation.
%
% Inputs:
%       p   = Hyperparameters, [a g1 g2 ... gd v], D=diag{g1,g2,...,gd}
%       M   = Input explanatory variable, of dimension {N0 x D0}
%       y   = Outcome (vector) from every input of the explantory variable
%       T   = Input explantory variable to be predicted for its outcome
%
% Outputs:
%       PR  = Prediction
%       SD  = Covariance matrix
%
%
% (C) Gaussian Process Toolbox 1.0
% (C) Copyright 2005-2007, Keith Neo <keith.neo@nuim.ie>
%                          http://www.hamilton.ie

[N0,D0]=size(M);
a=p(1);
g=p(2:D0+1);
v=p(D0+2);
Nt=size(T,1);
C=diag_add(eecov(M,M,a,g,N0),a*v,N0);
invC=inv(C);
Cy1=eecov(M,T,a,g,N0);
Ct1=eecov(T,T,a,g,Nt);

pr=Cy1'*invC*y;
sdev=Ct1-Cy1'*invC*Cy1;

function Q=eecov(M,N,a,g,N0)
Nt=size(N,1);
Z2=(repmat(permute(M,[1,3,2]),[1,Nt,1])-repmat(permute(N',[3,2 1]),[N0,1,1])).^2;
Q=a*exp(-0.5*sum(Z2.*repmat(permute(g',[3,2,1]),[N0,Nt,1]),3));

function Q=diag_add(Q,v,N0)
for i=1:N0
    Q(i,i)=Q(i,i)+v;
end