lrm_d_like.m

轨迹聚类工具箱及其相关说明和文章（台风路径分类等）。如有问题可与[e-mail]联系

function [Lhood,other] = lrm_d_like(varargin)
%LRM_D_LIKE  Calculate log-likelihood with LRM_D model.
%
%   [Lhood,Other] = LRM_D_LIKE(M,Trajs,[Options])
%    - M       : trained model
%    - Trajs   : 'Trajs' structure; See also CCTOOLBOX
%    - Options : see MODEL_LIKE
%
%   [Lhood,Other] = LRM_D_LIKE(M,X,Y,Seq,[Options])
%    - M       : trained model
%    - X,Y,Seq : curves in Sequence format; see CCTOOLBOX
%    - Options : see MODEL_LIKE
%
%   Other
%     .C   : classification labels

% Scott Gaffney   5 October 2003
% DataLab@UCI
% Department of Information and Computer Science
% University of California, Irvine, USA.

PROGNAME = 'lrm_d_like';
if (~nargin)
  try; help(PROGNAME); catch; end
  return;
end


%%% Handle Argument Processing
%%%
args = varargin; clear varargin;
n = length(args);
trajs=[]; X=[]; Y=[]; Seq=[]; Ops=[];
%
% Check for calling convention
%
% LRM_D_LIKE(M,Trajs,[Options])
if (n<4)
  M = args{1};
  trajs = args{2};
  if (n>2)
    Ops = args{3};
  end
  
% LRM_D_LIKE(M,X,Y,Seq,[Options])
else
  M = args{1};
  X = args{2};
  Y = args{3};
  Seq = args{4};
  if (n>4)
    Ops = args{5};
  end
end
%%
%%% End Argument Processing

if (isempty(Y))
  [Y,X,Seq] = trajs2seq(trajs,M.zero,M.Options.MinLen);
end
if (size(X,2)~=size(M.Mu,1))
  X = regmat(X,M.order);  
end


[N,D] = size(Y);
NumPoints = N*D;
K = M.K;
n = length(Seq)-1;
M.Mu = permute(M.Mu,[1 3 2]);

Pikd = zeros(n,K,D);
for k=1:K
  for d=1:D
    Mu       = M.Mu(:,d,k);
    sigma    = M.Sigma(k,d);
    s        = M.S(k,d);
    for i=1:n
      indx   = Seq(i):Seq(i+1)-1;
      ni    = length(indx);
      XMu    = X(indx,:)*Mu;
      iS     = eye(ni)/sigma - 1/(ni*sigma + sigma^2/s);
      Pikd(i,k,d) = mvnormpdf_inv(Y(indx,d)',XMu',iS);
    end
  end
end

scale = mean(mean(mean(Pikd)));
Pikd(:,:,1) = prod(Pikd./scale,3);
Pikd(:,:,1) = Pikd(:,:,1) .* (ones(n,1)*M.Alpha');
s = sum(Pikd(:,:,1),2);
z = find(s==0);
if (~isempty(z))
  Pikd(z,:,1) = realmin*1e100*(ones(length(z),1)*M.Alpha');
  s(z) = sum(Pikd(z,:,1),2);
end
Lhood = sum(log(s)) + n*D*log(scale);
Lhood = Lhood./NumPoints;

% Classify the sequences
if (nargout>1)
  [trash, other.C] = max(Pikd(:,:,1),[],2);
end