Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/loglikelihood_ssm.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 % Compute log-likelihood for SSM objects
4 %
5 % M Nofrarias 15-06-09
6 %
7 % $Id: loglikelihood_ssm.m,v 1.3 2011/11/16 08:52:49 nikos Exp $
8 %
9 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
10 function [loglk snr]= loglikelihood_ssm(xn,in,out,noise,model,params,inNames,outNames, spl, Amats, Bmats, Cmats, Dmats)
11
12 loglk = 0;
13 snr = 0;
14 switch class(in)
15 case 'ao'
16 % parameters
17 fs = noise(1).fs;
18 N = length(noise(1).y);
19
20 if (numel(in) == 1 && numel(out) == 1)
21 xn = double(xn);
22
23 spl = plist('set', 'for bode', ...
24 'outputs', outNames, ...
25 'inputs', inNames, ...
26 'reorganize', false,...
27 'f', in(1).x);
28
29 eval = copy(model,1);
30
31 % set parameter values
32 eval.doSetParameters(params, xn);
33
34 % make numeric
35 eval.doSubsParameters(params, true);
36
37 % do bode
38 h11 = bode(eval, spl, 'internal');
39
40 f = in.x;
41 % spectra to variance
42 C11 = (N*fs/2)*noise(1).y;
43
44 % compute elements of inverse cross-spectrum matrix
45 InvS11 = 1./C11;
46
47 % compute log-likelihood terms first, all at once does not cancel the
48 % imag part when multiplying x.*conj(x)
49 v1v1 = conj(out(1).y - h11.y.*in(1).y).*(out(1).y - h11.y.*in(1).y);
50
51 tmplt = h11.*in(1).y;
52
53 %computing SNR
54 snrexp = utils.math.stnr(tmplt,0,out(1).y,0,InvS11,0,0,0);
55
56 snr = snr + 20*log10(snrexp);
57
58 log1exp = sum(InvS11.*v1v1);
59
60 loglk = loglk + log1exp;
61
62 elseif (numel(in) == 2 && numel(out) == 2)
63
64 f = in(1).x;
65
66 xn = double(xn);
67
68 spl = plist('set', 'for bode', ...
69 'outputs', outNames, ...
70 'inputs', inNames, ...
71 'reorganize', false,...
72 'f', in(1).x);
73
74 eval = copy(model,1);
75
76 % set parameter values
77 eval.doSetParameters(params, xn);
78
79 % make numeric
80 eval.doSubsParameters(params);
81
82 % do bode
83 h = bode(eval, spl, 'internal');
84 h11 = h(1);
85 h12 = h(2);
86 h21 = h(3);
87 h22 = h(4);
88
89 % spectra to variance
90 C11 = (N*fs/2)*noise(1).y;
91 C22 = (N*fs/2)*noise(2).y;
92 C12 = (N*fs/2)*noise(3).y;
93 C21 = (N*fs/2)*noise(4).y;
94
95 % compute elements of inverse cross-spectrum matrix
96 InvS11 = (C22./(C11.*C22 - C12.*C21));
97 InvS22 = (C11./(C11.*C22 - C12.*C21));
98 InvS12 = (C21./(C11.*C22 - C12.*C21));
99 InvS21 = (C12./(C11.*C22 - C12.*C21));
100
101 % compute log-likelihood terms first, all at once does not cancel the
102 % imag part when multiplying x.*conj(x)
103 v1v1 = conj(out(1).y - h11.y.*in(1).y - h12.y.*in(2).y).*(out(1).y - h11.y.*in(1).y - h12.y.*in(2).y);
104 v2v2 = conj(out(2).y - h21.y.*in(1).y - h22.y.*in(2).y).*(out(2).y - h21.y.*in(1).y - h22.y.*in(2).y);
105 v1v2 = conj(out(1).y - h11.y.*in(1).y - h12.y.*in(2).y).*(out(2).y - h21.y.*in(1).y - h22.y.*in(2).y);
106 v2v1 = conj(out(2).y - h21.y.*in(1).y - h22.y.*in(2).y).*(out(1).y - h11.y.*in(1).y - h12.y.*in(2).y);
107
108 tmplt1 = h11.*in(1).y + h12.*in(2).y;
109 tmplt2 = h21.*in(1).y + h22.*in(2).y;
110
111 %computing SNR
112 snrexp = utils.math.stnr(tmplt1,tmplt2,out(1).y,out(2).y,InvS11,InvS22,InvS12,InvS21);
113
114 snr = snr + 20*log10(snrexp);
115
116 log1exp = sum(InvS11.*v1v1 + InvS22.*v2v2 - InvS12.*v1v2 - InvS21.*v2v1);
117
118 loglk = loglk + log1exp;
119
120 else
121 error('This method is only implemented for 1 input / 1 output model or for 2 inputs / 2 outputs models');
122 end
123
124 case 'matrix'
125 % parameters
126
127 fs = in(1).objs(1).fs;
128
129 % num. experiments
130 nexp = numel(in);
131
132 noutChannels = numel(out(1).objs);
133
134 N = length(noise(1).objs(1).y);
135
136 loglk = 0;
137 for nnn = 1:nexp
138
139 if ((numel(in(1).objs) == 1) && numel(out(1).objs) == 1)
140
141 freqs = in(nnn).objs(1).data.getX;
142
143 xn = double(xn);
144
145 spl = plist('set', 'for bode', ...
146 'outputs', outNames, ...
147 'inputs', inNames, ...
148 'reorganize', false,...
149 'f', freqs);
150
151
152 eval = copy(model(nnn),1);
153
154 % set parameter values
155 eval.doSetParameters(params, xn);
156
157 % make numeric
158 eval.doSubsParameters(params, true);
159
160 % do bode
161 h(:,1) = bode(eval, spl, 'internal');
162
163 for j = 1:noutChannels^2
164 % spectra to variance
165 % (N*fs/2)* this multiplication is done now in mcmc
166 C(:,j) = noise(nnn).objs(j).data.getY;
167 end
168
169 % compute elements of inverse cross-spectrum matrix
170 InvS11 = 1./C(:,1);
171
172
173 % compute log-likelihood terms first, all at once does not cancel the
174 % imag part when multiplying x.*conj(x)
175 in1 = in(nnn).objs(1).data.getY;
176 out1 = out(nnn).objs(1).data.getY;
177
178 % matrix index convention: H(1,1)->h(1) H(2,1)->h(2) H(1,2)->h(3) H(2,2)->h(4)
179 v1v1 = conj(out1 - h.getObjectAtIndex(1,1).y.*in1).*(out1 - h.getObjectAtIndex(1,1).y.*in1);
180
181 tmplt = h.getObjectAtIndex(1,1).y.*in1;
182
183 %computing SNR
184 snrexp = utils.math.stnr(tmplt,0,out1,0,InvS11,0,0,0);
185
186 snr = snr + 20*log10(snrexp);
187
188 log1exp = sum(InvS11.*v1v1);
189
190 loglk = loglk + log1exp;
191
192 elseif ((numel(in(1).objs) == 2) && numel(out(1).objs) == 2)
193
194 freqs = in(nnn).objs(1).data.getX;
195
196 xn = double(xn);
197
198 spl.pset('f', freqs);
199
200 eval = model(nnn);
201 eval.setA(Amats);
202 eval.setB(Bmats);
203 eval.setC(Cmats);
204 eval.setD(Dmats);
205
206 % set parameter values
207 eval.doSetParameters(params, xn);
208
209 % make numeric
210 eval.doSubsParameters(params, true);
211
212 % do bode
213 [h1 h2 h3 h4] = bode(eval, spl);
214
215 for j = 1:noutChannels^2
216 % spectra to variance
217 % (N*fs/2)* this multiplication is done now in mcmc
218 C(:,j) = noise(nnn).objs(j).data.getY;
219 end
220
221 % compute elements of inverse cross-spectrum matrix
222 detm = (C(:,1).*C(:,4) - C(:,2).*C(:,3));
223 InvS11 = C(:,4)./detm; %1 4
224 InvS22 = C(:,1)./detm; %4 1
225 InvS12 = C(:,2)./detm; %2 2
226 InvS21 = C(:,3)./detm; %3 3
227
228 % compute log-likelihood terms first, all at once does not cancel the
229 % imag part when multiplying x.*conj(x)
230 in1 = in(nnn).objs(1).data.getY;
231 in2 = in(nnn).objs(2).data.getY;
232 out1 = out(nnn).objs(1).data.getY;
233 out2 = out(nnn).objs(2).data.getY;
234
235 % matrix index convention: H(1,1)->h(1) H(2,1)->h(2) H(1,2)->h(3) H(2,2)->h(4)
236 v1v1 = conj(out1 - h1.*in1 - h3.*in2).*(out1 - h1.*in1 - h3.*in2);
237 v2v2 = conj(out2 - h2.*in1 - h4.*in2).*(out2 - h2.*in1 - h4.*in2);
238 v1v2 = conj(out1 - h1.*in1 - h3.*in2).*(out2 - h2.*in1 - h4.*in2);
239 v2v1 = conj(out2 - h2.*in1 - h4.*in2).*(out1 - h1.*in1 - h3.*in2);
240
241 tmplt1 = h1.*in1 + h3.*in2;
242 tmplt2 = h2.*in1 + h4.*in2;
243
244 %computing SNR
245 snrexp = utils.math.stnr(tmplt1,tmplt2,out1,out2,InvS11,InvS22,InvS12,InvS21);
246
247 snr = snr + 20*log10(snrexp);
248
249 log1exp = sum(InvS11.*v1v1 + InvS22.*v2v2 - InvS12.*v1v2 - InvS21.*v2v1);
250
251 loglk = loglk + log1exp;
252 else
253 error('This method is only implemented for 1 input / 1 output model or for 2 inputs / 2 outputs models');
254 end
255
256
257
258 end
259 end
260 end
261