view m-toolbox/classes/+utils/@math/loglikelihood.m @ 39:11e3ed9d2115
database-connection-manager
Implement databases listing in database connection dialog
author |
Daniele Nicolodi <nicolodi@science.unitn.it> |
date |
Mon, 05 Dec 2011 16:20:06 +0100 (2011-12-05) |
parents |
f0afece42f48 |
children |
|
line source
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute log-likelihood
%
% M Nofrarias 15-06-09
%
% $Id: loglikelihood.m,v 1.2 2011/03/24 19:59:50 ingo Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function loglk = loglikelihood(xn,in,out,noise,model,params)
% parameters
fs = noise(1).fs;
N = length(noise(1).y);
if numel(in) == 1 % 1 input / 1 output case
f = in(1).x;
% evaluate models
eval = model(1).setParams(params,xn);
h11 = double(eval);
% spectra to variance
C11 = (N*fs/2)*noise(1).y;
% compute elements of inverse cross-spectrum matrix
InvS11 = 1./C11;
% compute log-likelihood terms first, all at once does not cancel the
% imag part when multiplying x.*conj(x)
v1v1 = conj(out(1).y - h11.*in(1).y).*(out(1).y - h11.*in(1).y);
loglk = sum(InvS11.*v1v1);
elseif numel(in) == 2
if(numel(out) == 1) % 2 input / 1 output case
f = in(1).x;
% evaluate models
h11 = model(1).setParams(params,xn).double;
h12 = model(2).setParams(params,xn).double;
% spectra to variance
C11 = (N*fs/2)*noise(1).y;
% compute elements of inverse cross-spectrum matrix
InvS11 = 1./C11;
% compute log-likelihood terms first, all at once does not cancel the
% imag part when multiplying x.*conj(x)
in1 = in(1).data.getY;
f = in(1).data.getX;
in2 = in(2).data.getY;
out1 = out(1).data.getY;
v1v1 = conj(out1 - h11.*in1 - h12.*in2).*(out1 - h11.*in1 - h12.*in2);
loglk = sum(InvS11.*v1v1);
elseif(numel(out) == 2) % 2 input / 2 output case
f = in(1).x;
% evaluate models
h11 = model(1).setParams(params,xn).double;
h12 = model(2).setParams(params,xn).double;
h21 = model(3).setParams(params,xn).double;
h22 = model(4).setParams(params,xn).double;
% spectra to variance
C11 = (N*fs/2)*noise(1).y;
C22 = (N*fs/2)*noise(2).y;
C12 = (N*fs/2)*noise(3).y;
C21 = (N*fs/2)*noise(4).y;
% compute elements of inverse cross-spectrum matrix
InvS11 = (C22./(C11.*C22 - C12.*C21));
InvS22 = (C11./(C11.*C22 - C12.*C21));
InvS12 = (C21./(C11.*C22 - C12.*C21));
InvS21 = (C12./(C11.*C22 - C12.*C21));
% compute log-likelihood terms first, all at once does not cancel the
% imag part when multiplying x.*conj(x)
in1 = in(1).data.getY;
in2 = in(2).data.getY;
out1 = out(1).data.getY;
out2 = out(2).data.getY;
v1v1 = conj(out1 - h11.*in1 - h12.*in2).*(out1 - h11.*in1 - h12.*in2);
v2v2 = conj(out2 - h21.*in1 - h22.*in2).*(out2 - h21.*in1 - h22.*in2);
v1v2 = conj(out1 - h11.*in1 - h12.*in2).*(out2 - h21.*in1 - h22.*in2);
v2v1 = conj(out2 - h21.*in1 - h22.*in2).*(out1 - h11.*in1 - h12.*in2);
loglk = sum(InvS11.*v1v1 + InvS22.*v2v2 - InvS12.*v1v2 - InvS21.*v2v1);
end
end
end