www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@xregcovariance/mlelincost.m

    function [f,Beta,G,W,Li,T2] = mlelincost(cparams,c,y,X,W0,isNested,varargin)
%MLELINCOST MLE cost function for global covariance matrix with linear models
%
% [f,Beta,G,W,Li,T2] = MLELINCOST(cparams,c,y,X,W0,0)
% [f,Beta,G,W,Li,T2] = MLELINCOST(cparams,c,y,yhat,W0,1)

%  Copyright 2000-2006 The MathWorks, Inc. and Ford Global Technologies, Inc.

N = size(X,1);

% covariance model for Gamma
c = update(c,cparams);

% eval covariance function Wci*Wci'= inv(cov);
Wci = choltinv(c,W0);

if nargin<6 || isNested
    % Calculate weighted residuals Wc\(y-Zb);
    X = (Wci*X);
    y = Wci*y;
    if nargout==1 && N>1000 && issparse(X)
        % use q-less sparse qr decomp for speed
        [qy,R] = qr(X,y,0);
        % scalar result so divide by N
        r = sqrt((sum(y.*y)- sum(qy.*qy))/N);
    else
        Beta = X\y;
        r = y-X*Beta;
    end
else
    r = Wci*(y-X);
    Beta = [];
end

% log likelihood function
% log(det(W)) + (y-Zb)'/W*(y-Zb)
% e= 2*log(diag(Wc)) + r.^2;
if isempty(Wci)
    % diag is not well defined for empty matrices.
    wd = zeros(0,1);
else
    wd = diag(Wci);
end

f = sum( -2*log(wd) + r.^2);
% Note det(W)  = det(Wc)^2
%      det(Wc) = prod(diag(Wc)) as Wc is triangular
%      use log(prod(x)) = sum(log(x)) to avoid overflow

if nargout>1
    % don't calculate in optimisation

    % single unstructured covariance
    G = unstruct(c);

    Ng = size(G,1);

    W = cov(c,W0);

    % individual T^2 contributions
    T2 = reshape(r,Ng,length(r)/Ng);
    T2 = sum(T2.*T2,2);

    % individual likelihood contributions
    Li = reshape(-2*log(wd) + r.^2,Ng,length(r)/Ng);
    % sum over response features
    Li = sum(Li)';
end

drawnow