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

% GP2_IDEN_DERV         Calculates the derivative prediction and standard
% deviation of the posterior joint probability distribution for a Gaussian
% process.
%
% [PR,SD] = GP2_IDEN_DERV(p,M,y,T)
%
% This function computes the prediction and standard deviation of the
% derivative observation of the 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 of the derivative
%       SD  = Standard deviation for the derivative
%
%
% (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);

pr=zeros(Nt,D0);
sdev=zeros(Nt,D0);
for k=1:D0
    Z2=repmat(permute(T(:,k)',[3,2 1]),[N0,1,1])-repmat(permute(M(:,k),[1,3,2]),[1,Nt,1]);
    Cy11=-g(k)*Z2.*Cy1;
    pr(:,k)=Cy11'*(invC*y);
    sdev(:,k)=sqrt(a*g(k)-diag(Cy11'*invC*Cy11));
end

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