Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/loglikelihood.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % Compute log-likelihood
5 %
6 % M Nofrarias 15-06-09
7 %
8 % $Id: loglikelihood.m,v 1.2 2011/03/24 19:59:50 ingo Exp $
9 %
10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
11 function loglk = loglikelihood(xn,in,out,noise,model,params)
12
13 % parameters
14 fs = noise(1).fs;
15 N = length(noise(1).y);
16
17 if numel(in) == 1 % 1 input / 1 output case
18 f = in(1).x;
19 % evaluate models
20 eval = model(1).setParams(params,xn);
21 h11 = double(eval);
22
23 % spectra to variance
24 C11 = (N*fs/2)*noise(1).y;
25
26 % compute elements of inverse cross-spectrum matrix
27 InvS11 = 1./C11;
28
29 % compute log-likelihood terms first, all at once does not cancel the
30 % imag part when multiplying x.*conj(x)
31 v1v1 = conj(out(1).y - h11.*in(1).y).*(out(1).y - h11.*in(1).y);
32
33 loglk = sum(InvS11.*v1v1);
34
35 elseif numel(in) == 2
36
37 if(numel(out) == 1) % 2 input / 1 output case
38 f = in(1).x;
39
40 % evaluate models
41 h11 = model(1).setParams(params,xn).double;
42 h12 = model(2).setParams(params,xn).double;
43
44 % spectra to variance
45 C11 = (N*fs/2)*noise(1).y;
46
47 % compute elements of inverse cross-spectrum matrix
48 InvS11 = 1./C11;
49
50 % compute log-likelihood terms first, all at once does not cancel the
51 % imag part when multiplying x.*conj(x)
52 in1 = in(1).data.getY;
53 f = in(1).data.getX;
54 in2 = in(2).data.getY;
55 out1 = out(1).data.getY;
56
57 v1v1 = conj(out1 - h11.*in1 - h12.*in2).*(out1 - h11.*in1 - h12.*in2);
58
59 loglk = sum(InvS11.*v1v1);
60
61 elseif(numel(out) == 2) % 2 input / 2 output case
62 f = in(1).x;
63
64 % evaluate models
65 h11 = model(1).setParams(params,xn).double;
66 h12 = model(2).setParams(params,xn).double;
67 h21 = model(3).setParams(params,xn).double;
68 h22 = model(4).setParams(params,xn).double;
69
70 % spectra to variance
71 C11 = (N*fs/2)*noise(1).y;
72 C22 = (N*fs/2)*noise(2).y;
73 C12 = (N*fs/2)*noise(3).y;
74 C21 = (N*fs/2)*noise(4).y;
75
76 % compute elements of inverse cross-spectrum matrix
77 InvS11 = (C22./(C11.*C22 - C12.*C21));
78 InvS22 = (C11./(C11.*C22 - C12.*C21));
79 InvS12 = (C21./(C11.*C22 - C12.*C21));
80 InvS21 = (C12./(C11.*C22 - C12.*C21));
81
82 % compute log-likelihood terms first, all at once does not cancel the
83 % imag part when multiplying x.*conj(x)
84 in1 = in(1).data.getY;
85 in2 = in(2).data.getY;
86 out1 = out(1).data.getY;
87 out2 = out(2).data.getY;
88
89 v1v1 = conj(out1 - h11.*in1 - h12.*in2).*(out1 - h11.*in1 - h12.*in2);
90 v2v2 = conj(out2 - h21.*in1 - h22.*in2).*(out2 - h21.*in1 - h22.*in2);
91 v1v2 = conj(out1 - h11.*in1 - h12.*in2).*(out2 - h21.*in1 - h22.*in2);
92 v2v1 = conj(out2 - h21.*in1 - h22.*in2).*(out1 - h11.*in1 - h12.*in2);
93
94 loglk = sum(InvS11.*v1v1 + InvS22.*v2v2 - InvS12.*v1v2 - InvS21.*v2v1);
95 end
96 end
97
98 end
99