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

% GP2_IDEN_DURLEV       Calculates the prediction and standard deviation of
% the posterior joint probability distribution for the Gaussian process.
% The covariance matrices are strictly Toeplitz.
%
% [PR,SD] = GP2_IDEN_DURLEV(p,M,y,T)
%
% This function computes the prediction and the standard deviation 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.
%
% 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  = Standard deviation
%
% This function script utilises modified Durbin-Levinson's algorithm to
% construct fast algorithm for computation of task impossible by standard
% MATLAB command.
%
%
% (C) Gaussian Process Toeplitz Toolbox 1.0
% (C) Copyright 2005-2007, Keith Neo <keith.neo@nuim.ie>
%                          http://www.hamilton.ie

N0=size(M,1);
Nt=size(T,1);
a=p(1);
g=p(2);
v=p(3);
ct=toep2_genct(M,[1 g v]);
invb=toep2_invcy(ct,y);
if N0==Nt
    if sum(double(M==T))==N0
        pr=y-v*invb;
        ct=toep3_invdiag(ct);
        k=a*v;
        sdev=sqrt(k-k*v*ct);
    end
elseif Nt<2000
    C=diag_add(eecov_1d(M,M,a,g,N0),a*v,N0);
    invC=inv(C);
    Cy1=eecov_1d(M,T,a,g,N0);
    pr=(Cy1'*invb)/a;
    sdev=sqrt(a-diag(Cy1'*invC*Cy1));
else
    disp('Data size too big. Standard deviation result not computed');
    pr=zeros(Nt,1);
    for k=1:Nt
        tmp=exp(-0.5*g*(M'-T(k)).^2);
        pr(k)=tmp*invb;
    end
    sdev=0;
end

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

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