Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/loglikelihood_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 assuming a multivariate gaussian
4 % distribution
5 %
6 % INPUT
7 %
8 % - res, a vector of AOs containing residuals
9 % - noise, a vector of AOs containing noise
10 % - cutbefore, followed by the data samples to cut at the starting of the
11 % data series
12 % - cutafter, followed by the data samples to cut at the ending of the
13 % data series
14 %
15 % L Ferraioli 10-10-2010
16 %
17 % $Id: loglikelihood_td.m,v 1.1 2011/03/15 16:19:20 miquel Exp $
18 %
19 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
20 function loglk = loglikelihood_td(res,noise,varargin)
21
22 cutbefore = [];
23 cutafter = [];
24 if ~isempty(varargin)
25 for j=1:length(varargin)
26 if strcmp(varargin{j},'cutbefore')
27 cutbefore = varargin{j+1};
28 end
29 if strcmp(varargin{j},'cutafter')
30 cutafter = varargin{j+1};
31 end
32 end
33 end
34
35 nres = numel(res);
36 loglk = 0;
37
38 for ii=1:nres
39
40 yres = res(ii).y;
41 ynoise = noise(ii).y;
42 % willing to work with rows
43 if size(yres,2)<size(yres,1)
44 yres = yres.';
45 end
46 if size(ynoise,2)<size(ynoise,1)
47 ynoise = ynoise.';
48 end
49 if ~isempty(cutbefore)
50 yres(1:cutbefore) = [];
51 ynoise(1:cutbefore) = [];
52 end
53 if ~isempty(cutafter)
54 yres(end-cutafter:end) = [];
55 ynoise(end-cutafter:end) = [];
56 end
57
58 % R = utils.math.Rcovmat(yres);
59 %
60 % xx = R\yres.';
61 % Ndim = numel(yres);
62 %
63 % loglk = loglk + abs(xx'*xx)./Ndim;
64
65 cmat = utils.math.xCovmat(ynoise);
66
67 [L,p] = chol(cmat,'lower');
68 if p==0
69 xx = L\yres.';
70 Ndim = numel(yres);
71 else
72 q = p-1;
73 yyres = yres(1:q);
74
75 xx = L\yyres.';
76 Ndim = numel(yyres);
77 end
78
79 loglk = loglk + abs(xx'*xx)./Ndim;
80
81
82 end
83
84
85 %%% get cros terms
86 if nres>1
87
88 for ii=1:nres-1
89 for jj=ii+1:nres
90
91 yres1 = res(ii).y;
92 yres2 = res(jj).y;
93 ynoise1 = noise(ii).y;
94 ynoise2 = noise(jj).y;
95 % willing to work with rows
96 if size(yres1,2)<size(yres1,1)
97 yres1 = yres1.';
98 end
99 if size(yres2,2)<size(yres2,1)
100 yres2 = yres2.';
101 end
102 if size(ynoise1,2)<size(ynoise1,1)
103 ynoise1 = ynoise1.';
104 end
105 if size(ynoise2,2)<size(ynoise2,1)
106 ynoise2 = ynoise2.';
107 end
108 if ~isempty(cutbefore)
109 yres1(1:cutbefore) = [];
110 yres2(1:cutbefore) = [];
111 ynoise1(1:cutbefore) = [];
112 ynoise2(1:cutbefore) = [];
113 end
114 if ~isempty(cutafter)
115 yres1(end-cutafter:end) = [];
116 yres2(end-cutafter:end) = [];
117 ynoise1(end-cutafter:end) = [];
118 ynoise2(end-cutafter:end) = [];
119 end
120
121 % Rx = utils.math.Rcovmat(yres1);
122 % Ry = utils.math.Rcovmat(yres2);
123 %
124 % xx = Rx\yres1.';
125 % yy = Ry\yres2.';
126 %
127 % Ndim = numel(yres1);
128 %
129 % loglk = loglk + 2.*abs(xx'*yy)./Ndim;
130
131 cmat = utils.math.xCovmat(ynoise1,ynoise2);
132
133 [L,p] = chol(cmat,'lower');
134 if p==0
135 xx = L\yres1.';
136 Ndim = numel(yres1);
137 yy = L\yres2.';
138 else
139 q = p-1;
140 yyres1 = yres1(1:q);
141 yyres2 = yres2(1:q);
142
143 xx = L\yyres1.';
144 Ndim = numel(yyres1);
145 yy = L\yyres2.';
146 end
147
148 loglk = loglk + 2.*abs(xx'*yy)./Ndim;
149
150
151 end
152 end
153
154 end
155
156
157 end