www.pudn.com > WebsiteCode.zip > gelikelihood.m, change:2008-07-24,size:3264b

```function [NegLnLike,PsiDot,U]=gelikelihood(x)
%[NegLnLike,PsiDot,U]=gelikelihood(x,ModelInfo)
%
% Calculates the negative of the concentrated log-likelihood
% Inputs:
%	x - vetor of paramters, either 1 x k specifies theta only or
%	1 x k+1 specifies theta and lambda
%	ModelInfo.X - n x n matrix of sample locations
%	ModelInfo.y - 1 x (2k+1)n vector of observed data
%
% Outputs:
%	NegLnLike - concentrated log-likelihood *-1 for minimising
%	PsiDot - differentiated correlation matrix
%	L,U - LU decomposition of PsiDot
%
% Copyright 2007 A I J Forrester
%
% This program is free software: you can redistribute it and/or modify  it
% the Free Software Foundation, either version 3 of the License, or 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 Lesser
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License and GNU
% Lesser General Public License along with this program. If not, see

global ModelInfo
X=ModelInfo.X;
y=ModelInfo.y;
theta=10.^x;
n=size(X,1);
k=size(X,2);
one=[ones(n,1); zeros(n*k,1)];

% build upper half of correlation matrix
% pre-allocate memory
Psi=zeros(n,n);
for i=1:n
for j=i+1:n
Psi(i,j)=exp(-sum(theta.*(X(i,:)-X(j,:)).^2)); % build top right corner
end
end
Psi=Psi+Psi'+eye(n)+eye(n).*eps;
% build differentiated matrices
% Pre-allocate memory
PsiDot=zeros((k+1)*n,(k+1)*n);
% Build upper half of PsiDot
for l=1:k
for m=l:k
if l==1
% Build upper half of dPsidX
for i=1:n
for j=i+1:n
PsiDot(i,m*n+j)=2*theta(m)*(X(i,m)-X(j,m))*Psi(i,j);
end
end
% Add upper and lower halves
PsiDot(1:n,m*n+1:(m+1)*n)=PsiDot(1:n,m*n+1:(m+1)*n)...
-PsiDot(1:n,m*n+1:(m+1)*n)';
end
if m==l
% Build upper half of d2PsidX^2
for i=1:n
for j=i+1:n
PsiDot(l*n+i,m*n+j)=...
(2*theta(l)-4*theta(l)^2*(X(i,l)-X(j,l))^2)*Psi(i,j);
end
end
PsiDot(l*n+1:(l+1)*n,m*n+1:(m+1)*n)=PsiDot(l*n+1:(l+1)*n,m*n+1:(m+1)*n)...
+eye(n).*theta(l);
else
% Build upper half of d2PsidXdX
for i=1:n
for j=i+1:n
PsiDot(l*n+i,m*n+j)=-4*theta(l)*theta(m)*(X(i,l)-X(j,l))*(X(i,m)-X(j,m))*Psi(i,j);
end
end
% Add upper and lower halves
PsiDot(l*n+1:(l+1)*n,m*n+1:(m+1)*n)=PsiDot(l*n+1:(l+1)*n,m*n+1:(m+1)*n)...
+PsiDot(l*n+1:(l+1)*n,m*n+1:(m+1)*n)';
end
end
end
% Add upper and lower halves to Psi
PsiDot=[Psi zeros(n,k*n);zeros(k*n,(k+1)*n)]+PsiDot+PsiDot';

%%% concentrated log-likelihood calculation %%%
% Cholesky factorisation
[U,p]=chol(PsiDot);
if p>0
NegLnLike=1e4;
else

% twice sum logs of diagonal to find ln(abs(det(Psi)))
lndetPsi=2*sum(log(abs(diag(U))));
%%% use back-substitution of LU decomposition instead of inverse
mu=(one'*(U\(U'\y)))/(one'*(U\(U'\one)));
sigma_sqr=((y-one*mu)'*(U\(U'\(y-one*mu))))/((k+1)*n);
NegLnLike=-1*(-(((k+1)*n)/2)*log(sigma_sqr) - 0.5*lndetPsi);
end
```