Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/chisquare_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 %
3 % Compute log-likelihood in time domain for SSM objects
4 %
5 % INPUT
6 %
7 % - in, a vector of input signals aos
8 % - out, a vector of output data aos
9 % - parvals, a vector with parameters values
10 % - parnames, a cell array with parameters names
11 % - model, an ssm model
12 % - inNames, A cell-array of input port names corresponding to the
13 % different input AOs
14 % - outNames, A cell-array of output ports to return
15 % - cutbefore, followed by the data samples to cut at the starting of the
16 % data series
17 % - cutafter, followed by the data samples to cut at the ending of the
18 % data series
19 %
20 % L Ferraioli 10-10-2010
21 %
22 % $Id: chisquare_ssm_td.m,v 1.1 2010/11/16 16:41:37 luigi Exp $
23 %
24 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
25 function chi2 = chisquare_ssm_td(xp,in,out,parnames,model,inNames,outNames,varargin)
26 % xn,in,out,noise,model,params,inNames,outNames
27 cutbefore = [];
28 cutafter = [];
29 if ~isempty(varargin)
30 for j=1:length(varargin)
31 if strcmp(varargin{j},'cutbefore')
32 cutbefore = varargin{j+1};
33 end
34 if strcmp(varargin{j},'cutafter')
35 cutafter = varargin{j+1};
36 end
37 end
38 end
39
40 xp = double(xp);
41 fs = out(1).fs;
42
43 % set parameters in the model
44 evalm = model.setParameters(plist('names',parnames,'values',xp));
45 evalm.keepParameters();
46 evalm.modifyTimeStep(plist('newtimestep',1/fs));
47
48 %%% get expected outputs
49 plsym = plist('AOS VARIABLE NAMES',inNames,...
50 'RETURN OUTPUTS',outNames,...
51 'AOS',in);
52 eo = simulate(evalm,plsym);
53
54 %%% get measurement noise
55 res = out-eo;
56 resy = res.y;
57 if size(resy,2)>size(resy,1)
58 resy = resy.';
59 end
60
61 if ~isempty(cutbefore)
62 resy(1:cutbefore) = [];
63 end
64 if ~isempty(cutafter)
65 resy(end-cutafter:end) = [];
66 end
67
68 chi2 = resy'*resy;
69 chi2 = chi2/(numel(resy)-numel(xp));
70
71
72
73
74 end