Mercurial > hg > ltpda

comparison m-toolbox/test/utils/test_loglikehood_ssm_td.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 % Test loglikehood_ssm_td
3 %
4 % Test with an harmonic hoscillator and simplex parameters search
5 %
6 % L Ferraioli 11-10-2010
7 %
8 % $Id: test_loglikehood_ssm_td.m,v 1.3 2010/11/26 13:52:32 luigi Exp $
9 %
10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
11 %% set parameters
12
13 m = 2; % kg
14 k = 0.1; % kg s^-2
15 damp = 0.05; % kg s^-1
16
17 fAmp = 3;
18 NAmp = 1;
19 fs = 1;
20
21 cutbefore = 100*fs;
22 cutafter = 500*fs;
23
24 f = logspace(-5+log10(fs),log10(fs/2),300);
25
26 %% Load model
27
28 sys = ssm(plist('built-in', 'HARMONIC_OSC_1D', 'SYMBOLIC PARAMS',{'M', 'K', 'VBETA'}));
29 sys = sys.setParameters(plist('names',{'M', 'K', 'VBETA'},'values',[m k damp]));
30 sys.keepParameters();
31 sys.modifyTimeStep(plist('newtimestep',1/fs));
32
33 pl = plist('inputs', 'COMMAND.Force', 'HARMONIC_OSC_1D.outputs', 'Position', 'f', f);
34 h = bode(sys,pl);
35
36 iplot(h)
37
38 %% input force and readout noise
39
40 % input force
41 fc = ao(plist('waveform','sine wave','A',fAmp,'f',1e-1,'phi',0,'toff',cutbefore/fs,'nsecs',6e2,'fs',fs));
42 in = zeropad(fc,plist('Position','post','N',cutafter));
43
44
45 rn = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', in.nsecs));
46 rn = rn.*NAmp;
47
48
49
50 %% Simulate the system
51
52
53 plsym = plist('AOS VARIABLE NAMES',{'COMMAND.Force','noise.readout'},...
54 'RETURN OUTPUTS','HARMONIC_OSC_1D.position',...
55 'AOS',[in,rn]);
56 out = simulate(sys,plsym);
57
58 iplot(out,rn)
59
60
61
62
63 %% Do chisquare parameter search
64
65
66 mod = ssm(plist('built-in', 'HARMONIC_OSC_1D',...
67 'SYMBOLIC PARAMS',{'M', 'K', 'VBETA'}));
68 inNames = {'COMMAND.Force'};
69 outNames = {'HARMONIC_OSC_1D.position'};
70 parnames = {'m', 'k', 'vbeta'};
71
72
73 x0 = [1 1 1];
74 x0 = [2.5 0.08 0.1];
75
76
77 % loglk =
78 % loglikehood_ssm_td(xp,in,out,parnames,model,inNames,outNames,varargin)
79
80 xp2 = fminsearch(@(x)utils.math.chisquare_ssm_td(x,in,out,parnames,mod,inNames,outNames,'cutbefore',cutbefore,'cutafter',cutafter),x0,optimset('Display','iter'));
81
82
83
84 %%
85
86 sys = copy(mod,1);
87 sys = sys.setParameters(plist('names',{'M', 'K', 'VBETA'},'values',xp2));
88 sys.keepParameters();
89 sys.modifyTimeStep(plist('newtimestep',1/fs));
90
91 plsym = plist('AOS VARIABLE NAMES',{'COMMAND.Force'},...
92 'RETURN OUTPUTS','HARMONIC_OSC_1D.position',...
93 'AOS',[in]);
94 outfit2 = simulate(sys,plsym);
95
96 iplot(out,outfit2)
97 iplot(out-outfit2)
98 iplot(out-outfit2-rn)
99
100 % xp =
101 %
102 % 2.22352087931512 0.588010809495998 0.195689637008003
103
104
105
106 %% Do parameters search
107
108
109 mod = ssm(plist('built-in', 'HARMONIC_OSC_1D',...
110 'SYMBOLIC PARAMS',{'M', 'K', 'VBETA'}));
111 inNames = {'COMMAND.Force'};
112 outNames = {'HARMONIC_OSC_1D.position'};
113 inNoiseNames = {'noise.readout'};
114 inNoise = rn;
115 outNoiseNames = {'HARMONIC_OSC_1D.position'};
116 parnames = {'M', 'K', 'VBETA'};
117
118
119 % x0 = xp2;
120 x0 = [2.5 0.08 0.1];
121
122 % loglk =
123 % loglikehood_ssm_td(xp,in,out,parnames,model,inNames,outNames,varargin)
124
125
126
127 xp = fminsearch(@(x)utils.math.loglikehood_ssm_td(x,in,out,parnames,mod,inNames,outNames,rn,'cutbefore',cutbefore,'cutafter',cutafter),x0,optimset('Display','iter'));
128
129
130
131 %%
132
133 sys = copy(mod,1);
134 sys = sys.setParameters(plist('names',{'M', 'K', 'VBETA'},'values',xp));
135 sys.keepParameters();
136 sys.modifyTimeStep(plist('newtimestep',1/fs));
137
138 plsym = plist('AOS VARIABLE NAMES',{'COMMAND.Force'},...
139 'RETURN OUTPUTS','HARMONIC_OSC_1D.position',...
140 'AOS',[in]);
141 outfit = simulate(sys,plsym);
142
143 iplot(out,outfit)
144 iplot(out-outfit)
145 iplot(out-outfit-rn)
146
147
148 %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
149 %% test covariace estimation
150
151 mod = ssm(plist('built-in', 'HARMONIC_OSC_1D',...
152 'SYMBOLIC PARAMS',{'M', 'K', 'VBETA'}));
153 inNames = {'COMMAND.Force'};
154 outNames = {'HARMONIC_OSC_1D.position'};
155 parnames = {'m', 'k', 'vbeta'};
156
157 sys = copy(mod,1);
158 sys = sys.setParameters(plist('names',{'M', 'K', 'VBETA'},'values',[m k damp]));
159 sys.keepParameters();
160 sys.modifyTimeStep(plist('newtimestep',1/fs));
161
162 plsym = plist('AOS VARIABLE NAMES',{'COMMAND.Force'},...
163 'RETURN OUTPUTS','HARMONIC_OSC_1D.position',...
164 'AOS',[in]);
165 modeout = simulate(sys,plsym);
166
167 % get noise
168 noise = out - modeout;
169 nn = noise.y;
170 % get covariance matrix
171 cmat = utils.math.xCovmat(nn);
172 cmatseq = cmat(1,:);
173
174 % get covariance sequance
175 [cseq,lags] = xcov(nn,'unbiased');
176 cseq = ifftshift(cseq);
177 lags=ifftshift(lags);
178 cseqr = cseq(1:numel(cmatseq));
179 lagsr = lags(1:numel(cmatseq));
180 lagsr = lagsr.';
181
182 % compare results
183 figure
184 plot(lagsr,cmatseq.','k')
185 hold on
186 plot(lagsr,cseqr,'r')
187
188 figure
189 plot(lagsr,cmatseq.'-cseqr,'k')
190
191 [R,p] = chol(cmat);
192 d = eig(cmat);
193 % md = min(d);
194 % cmat2 = cmat + abs(md).*eye(size(cmat));
195 % [R,p] = chol(cmat2);
196
197 ccmat = fliplr(flipud(cmat));
198 [RR,pp] = chol(ccmat);
199
200 %% correlation matrix
201
202 [corr,sig] = utils.math.cov2corr(cmat);
203
204 % test positive definiteness
205 [R,p] = chol(corr);
206
207 % get a valid correlation
208 newcorr = validcorr(corr);
209
210 %%
211
212 % get covariance sequance
213 [cseq,lags] = xcov(nn,'unbiased');
214 Num = numel(nn);
215 for ii=1:Num
216
217 ccmat(ii,:) = cseq(Num-ii+1:2*Num-ii);
218
219 end
220
221 [R,p] = chol(ccmat);
222
223
224 %% test the likelihood
225
226 mod = ssm(plist('built-in', 'HARMONIC_OSC_1D',...
227 'SYMBOLIC PARAMS',{'M', 'K', 'VBETA'}));
228 inNames = {'COMMAND.Force'};
229 outNames = {'HARMONIC_OSC_1D.position'};
230 parnames = {'m', 'k', 'vbeta'};
231
232 sys = copy(mod,1);
233 sys = sys.setParameters(plist('names',{'M', 'K', 'VBETA'},'values',[m k damp]));
234 sys.keepParameters();
235 sys.modifyTimeStep(plist('newtimestep',1/fs));
236
237 plsym = plist('AOS VARIABLE NAMES',{'COMMAND.Force'},...
238 'RETURN OUTPUTS','HARMONIC_OSC_1D.position',...
239 'AOS',[in]);
240 modeout1 = simulate(sys,plsym);
241
242 sys = copy(mod,1);
243 sys = sys.setParameters(plist('names',{'M', 'K', 'VBETA'},'values',[1 1 1]));
244 sys.keepParameters();
245 sys.modifyTimeStep(plist('newtimestep',1/fs));
246
247 plsym = plist('AOS VARIABLE NAMES',{'COMMAND.Force'},...
248 'RETURN OUTPUTS','HARMONIC_OSC_1D.position',...
249 'AOS',[in]);
250 modeout2 = simulate(sys,plsym);
251
252 % xp = [1.1 0.6 0.01];
253
254 sys = copy(mod,1);
255 sys = sys.setParameters(plist('names',{'M', 'K', 'VBETA'},'values',xp));
256 sys.keepParameters();
257 sys.modifyTimeStep(plist('newtimestep',1/fs));
258
259 plsym = plist('AOS VARIABLE NAMES',{'COMMAND.Force'},...
260 'RETURN OUTPUTS','HARMONIC_OSC_1D.position',...
261 'AOS',[in]);
262 modeout3 = simulate(sys,plsym);
263
264 % xp2 = [2.22352087931512 0.588010809495998 0.195689637008003];
265
266 sys = copy(mod,1);
267 sys = sys.setParameters(plist('names',{'M', 'K', 'VBETA'},'values',xp2));
268 sys.keepParameters();
269 sys.modifyTimeStep(plist('newtimestep',1/fs));
270
271 plsym = plist('AOS VARIABLE NAMES',{'COMMAND.Force'},...
272 'RETURN OUTPUTS','HARMONIC_OSC_1D.position',...
273 'AOS',[in]);
274 modeout4 = simulate(sys,plsym);
275
276 % get noise
277 noise1 = out - modeout1; % true values
278 noise2 = out - modeout2; % starting guess
279 noise3 = out - modeout3; % fit values
280 noise4 = out - modeout4; % chi2 fit values
281
282
283 lik1 = utils.math.loglikehood_td(noise1,'cutbefore',cutbefore,'cutafter',cutafter)
284 lik2 = utils.math.loglikehood_td(noise2,'cutbefore',cutbefore,'cutafter',cutafter)
285 lik3 = utils.math.loglikehood_td(noise3,'cutbefore',cutbefore,'cutafter',cutafter)
286 lik4 = utils.math.loglikehood_td(noise4,'cutbefore',cutbefore,'cutafter',cutafter)
287
288 %% test a uniform grid
289
290 mod = ssm(plist('built-in', 'HARMONIC_OSC_1D',...
291 'SYMBOLIC PARAMS',{'M', 'K', 'VBETA'}));
292 inNames = {'COMMAND.Force'};
293 outNames = {'HARMONIC_OSC_1D.position'};
294 parnames = {'m', 'k', 'vbeta'};
295
296 mm = [1.8:0.02:2.2];
297 kc = [0.1:0.02:1];
298 dd = [0.01:0.02:1];
299
300 idx = 1;
301 pars = struct;
302
303 Ntot = numel(mm)*numel(kc)*numel(dd);
304
305 lik = zeros(Ntot,1);
306 %%
307
308 for ii=1:numel(mm)
309 for jj=1:numel(kc)
310 for kk=1:numel(dd)
311
312 fprintf('step %s of %s\n',num2str(idx),num2str(Ntot));
313 pars(idx).m = mm(ii);
314 pars(idx).k = kc(jj);
315 pars(idx).damp = dd(kk);
316
317 sys = copy(mod,1);
318 sys = sys.setParameters(plist('names',{'M', 'K', 'VBETA'},'values',[mm(ii) kc(jj) dd(kk)]));
319 sys.keepParameters();
320 sys.modifyTimeStep(plist('newtimestep',1/fs));
321
322 plsym = plist('AOS VARIABLE NAMES',{'COMMAND.Force'},...
323 'RETURN OUTPUTS','HARMONIC_OSC_1D.position',...
324 'AOS',[in]);
325 modeout = simulate(sys,plsym);
326
327 noise = out - modeout;
328
329 lik(idx) = utils.math.loglikehood_td(noise,'cutbefore',cutbefore,'cutafter',cutafter);
330
331 idx = idx +1;
332 end
333 end
334 end
335