Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/confint.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 % CONFINT Calculates confidence levels and variance for psd, lpsd, cohere, lcohere and curvefit parameters
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: CONFINT Input psd, mscohere (magnitude square coherence)
5 % and return confidence levels and variance for them.
6 % Spectra are assumed to be calculated with the WOSA method (Welch's
7 % Overlapped Segment Averaging Method)
8 %
9 % CALL: out = confint(a,pl)
10 %
11 %
12 % INPUTS:
13 % a - input analysis objects containing power spectral
14 % densities or magintude squared coherence.
15 % pl - input parameter list
16 %
17 % OUTPUTS:
18 % out - a collection object containing:
19 % lcl - lower confidence level
20 % ucl - upper confidence level
21 % var - expected spectrum variance
22 %
23 %
24 % If the last input argument is a parameter list (plist).
25 % The following parameters are recognised.
26 %
27 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'confint')">Parameters Description</a>
28 %
29 %
30 %
31 %
32 %
33 % VERSION: $Id: confint.m,v 1.20 2011/04/29 13:54:36 luigi Exp $
34 %
35 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36
37 function varargout = confint(varargin)
38
39 %%% Check if this is a call for parameters
40 if utils.helper.isinfocall(varargin{:})
41 varargout{1} = getInfo(varargin{3});
42 return
43 end
44
45 import utils.const.*
46 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
47
48 %%% Collect input variable names
49 in_names = cell(size(varargin));
50 for ii = 1:nargin,in_names{ii} = inputname(ii);end
51
52 %%% Collect all AOs
53 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
54 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
55
56 %%% avoid multiple AO at input
57 if numel(as)>1
58 error('!!! Too many input AOs, CONFINT can process only one AO per time !!!')
59 end
60
61 %%% avoid input modification
62 if nargout == 0
63 error('!!! CONFINT cannot be used as a modifier. Please give an output variable !!!');
64 end
65
66 %%% Parse plists
67 pl = parse(pl, getDefaultPlist());
68
69 %%% Find parameters
70 mtd = lower(find(pl, 'method'));
71 conf = find(pl, 'conf');
72 dof = find(pl, 'dof');
73 conf = conf/100; % go from percentage to fractional
74 Ntot = find(pl,'DataLength');
75
76 %%% check that fsdata is input
77 if ~isa(as.data, 'fsdata')
78 error('!!! Non-fsdata input, CONFINT can process only fsdata !!!')
79 end
80
81
82 % looking to dof
83 if isempty(dof)
84 calcdof = true;
85 else
86 if isa(dof, 'ao')
87 dof = dof.data.y;
88 calcdof = false;
89 else
90 calcdof = false;
91 end
92 end
93
94
95 %%% switching over methods
96 switch mtd
97 case 'psd'
98 %%% confidence levels for spectra calculated with psd
99
100 % calculating dof
101 if calcdof
102 dofs = getdof(as,plist('method',mtd,'DataLength',Ntot));
103 dof = dofs.y;
104 end % if calcdof
105 dof = round(dof);
106 if length(dof)~=1
107 error('!!! CONFINT for ao/psd method, dof must be a single number')
108 end
109
110 % Calculating Confidence Levels factors
111 alfa = 1 - conf;
112 c = utils.math.Chi2inv([1-alfa/2 alfa/2],dof);
113 c = dof./c;
114
115 % calculating variance
116 expvar = ((as.data.y).^2).*2./dof;
117
118 % calculating confidence levels
119 lwb = as.data.y.*c(1);
120 upb = as.data.y.*c(2);
121
122 case 'lpsd'
123 %%% confidence levels for spectra calculated with lpsd
124
125 % calculating dof
126 if calcdof
127 dofs = getdof(as,plist('method',mtd,'DataLength',Ntot));
128 dofs = dofs.y;
129
130 % extract number of frequencies bins
131 nf = length(as.x);
132
133 cl = ones(nf,2);
134 for jj=1:nf
135
136 % Calculating Confidence Levels factors
137 alfa = 1 - conf;
138 c = utils.math.Chi2inv([1-alfa/2 alfa/2],dofs(jj));
139 c = dofs(jj)./c;
140
141 % storing c
142 cl(jj,1) = c(1);
143 cl(jj,2) = c(2);
144 end % for jj=1:nf
145 else % if calcdof
146 if length(dof)~=length(as.x)
147 error('!!! CONFINT for ao/lpsd method, dof must be a vector of the same length of the frequencies vector')
148 end
149 dofs = round(dof);
150 cl = ones(length(as.x),2);
151 for jj = 1:length(as.x)
152 % Calculating Confidence Levels factors
153 alfa = 1 - conf;
154 c = utils.math.Chi2inv([1-alfa/2 alfa/2],dofs(jj));
155 c = dofs(jj)./c;
156
157 % storing c
158 cl(jj,1) = c(1);
159 cl(jj,2) = c(2);
160 end
161 end % if calcdof
162 % willing to work with columns
163 dy = as.data.y;
164 [ii,kk] = size(dy);
165 if ii<kk
166 dy = dy.';
167 rsp = true;
168 else
169 rsp = false;
170 end
171 % calculating variance
172 expvar = ((dy).^2).*2./dofs;
173
174 % calculating confidence levels
175 lwb = dy.*cl(:,1);
176 upb = dy.*cl(:,2);
177
178 % reshaping if necessary
179 if rsp
180 expvar = expvar.';
181 lwb = lwb.';
182 upb = upb.';
183 end
184
185 case 'mscohere'
186 %%% confidence levels for mscohere calculated with ao/cohere
187
188 % calculating dof
189 if calcdof
190 dofs = getdof(as,plist('method',mtd,'DataLength',Ntot));
191 dof = dofs.y;
192 end % if calcdof
193 dof = round(dof);
194 if length(dof)~=1
195 error('!!! CONFINT for ao/cohere method, dof must be a single number')
196 end
197
198 % Defining Y variable
199 Y = atanh(sqrt(as.data.y));
200
201 % Calculating Confidence Levels factor
202 alfa = 1 - conf;
203 c = -sqrt(2).*erfcinv(2*(1-alfa/2))./sqrt(dof);
204 Ylwb = Y - c;
205 Yupb = Y + c;
206
207 % calculating confidence levels
208 lwb = tanh(Ylwb).^2;
209 upb = tanh(Yupb).^2;
210
211 % calculating variance
212 expvar = ((1-(as.data.y).^2).^2).*((as.data.y).^2).*4./dof;
213
214 case 'mslcohere'
215 %%% confidence levels for spectra calculated with lpsd
216
217 % calculating dof
218 if calcdof
219 dofs = getdof(as,plist('method',mtd,'DataLength',Ntot));
220 dofs = dofs.y;
221
222 % extract number of frequencies bins
223 nf = length(as.x);
224
225 % willing to work with columns
226 dy = as.data.y;
227 [ii,kk] = size(dy);
228 if ii<kk
229 dy = dy.';
230 rsp = true;
231 else
232 rsp = false;
233 end
234
235 % Defining Y variable
236 Y = atanh(sqrt(dy));
237
238 cl = ones(nf,2);
239 for jj=1:nf
240
241 % Calculating Confidence Levels factors
242 alfa = 1 - conf;
243 c = -sqrt(2).*erfcinv(2*(1-alfa/2))./sqrt(dofs(jj));
244
245 % storing c and dof
246 cl(jj,1) = Y(jj) - c;
247 cl(jj,2) = Y(jj) + c;
248 end % for jj=1:nf
249
250 else % if calcdof
251 if length(dof)~=length(as.x)
252 error('!!! CONFINT for ao/lcohere method, dof must be a vector of the same length of the frequencies vector')
253 end
254 dofs = round(dof);
255
256 % willing to work with columns
257 dy = as.data.y;
258 [ii,kk] = size(dy);
259 if ii<kk
260 dy = dy.';
261 rsp = true;
262 else
263 rsp = false;
264 end
265
266 % Defining Y variable
267 Y = atanh(sqrt(dy));
268
269 cl = ones(length(as.x),2);
270 for jj = 1:length(as.x)
271 % Calculating Confidence Levels factors
272 alfa = 1 - conf;
273 c = -sqrt(2).*erfcinv(2*(1-alfa/2))./sqrt(dofs(jj));
274
275 % storing c
276 cl(jj,1) = Y(jj) - c;
277 cl(jj,2) = Y(jj) + c;
278 end
279 end % if calcdof
280
281 % calculating variance
282 expvar = ((1-(dy).^2).^2).*((dy).^2).*4./dofs;
283
284 % get not well defined coherence estimations
285 idd = dofs<=2;
286
287 % calculating confidence levels
288 lwb = tanh(cl(:,1));
289 % remove negative elements
290 idx = lwb < 0;
291 lwb(idx) = 0;
292 % set lower bound to zero in points where coharence is not well defined
293 lwb(idd) = 0;
294 upb = tanh(cl(:,2));
295 % set upper bound to one in points where coharence is not well
296 % defined
297 upb(idd) = 1;
298 lwb = lwb.^2;
299 upb = upb.^2;
300
301 % reshaping if necessary
302 if rsp
303 expvar = expvar.';
304 lwb = lwb.';
305 upb = upb.';
306 end
307
308
309
310 end %switch mtd
311
312 % Output data
313
314
315 % defining units
316 inputunit = get(as.data,'yunits');
317 varunit = unit(inputunit.^2);
318 varunit.simplify;
319 levunit = inputunit;
320 levunit.simplify;
321
322
323 % variance
324 plvar = plist('xvals', as.data.x, 'yvals', expvar, 'type', 'fsdata');
325 ovar = ao(plvar);
326
327 ovar.setFs(as.data.fs);
328 ovar.setT0(as.data.t0);
329 ovar.data.setEnbw(as.data.enbw);
330 ovar.data.setNavs(as.data.navs);
331 ovar.setXunits(as.data.xunits);
332 ovar.setYunits(varunit);
333 % Set output AO name
334 ovar.name = sprintf('var(%s)', ao_invars{:});
335
336 % lower confidence level
337 pllwb = plist('xvals', as.data.x, 'yvals', lwb, 'type', 'fsdata');
338 olwb = ao(pllwb);
339
340 olwb.setFs(as.data.fs);
341 olwb.setT0(as.data.t0);
342 olwb.data.setEnbw(as.data.enbw);
343 olwb.data.setNavs(as.data.navs);
344 olwb.setXunits(copy(as.data.xunits,1));
345 olwb.setYunits(levunit);
346 % Set output AO name
347 clev = [num2str(conf*100) '%'];
348 olwb.name = sprintf('%s_low_conf_level(%s)', clev, ao_invars{:});
349
350 % upper confidence level
351 plupb = plist('xvals', as.data.x, 'yvals', upb, 'type', 'fsdata');
352 oupb = ao(plupb);
353
354 oupb.setFs(as.data.fs);
355 oupb.setT0(as.data.t0);
356 oupb.data.setEnbw(as.data.enbw);
357 oupb.data.setNavs(as.data.navs);
358 oupb.setXunits(copy(as.data.xunits,1));
359 oupb.setYunits(levunit);
360 % Set output AO name
361 oupb.name = sprintf('%s_up_conf_level(%s)', clev, ao_invars{:});
362
363
364 outobj = collection(olwb,oupb,ovar);
365 outobj.setName(sprintf('%s conf levels for %s', clev, ao_invars{:}));
366 outobj.addHistory(getInfo('None'), pl, [ao_invars(:)], [as.hist]);
367
368 varargout{1} = outobj;
369
370
371 end
372
373 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
374 % Local Functions %
375 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
376
377
378 %----- LTPDA FUNCTIONS ----------------------------------------------------
379 %--------------------------------------------------------------------------
380 % Get Info Object
381 %--------------------------------------------------------------------------
382 function ii = getInfo(varargin)
383 if nargin == 1 && strcmpi(varargin{1}, 'None')
384 sets = {};
385 pl = [];
386 else
387 sets = {'Default'};
388 pl = getDefaultPlist;
389 end
390 % Build info object
391 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: confint.m,v 1.20 2011/04/29 13:54:36 luigi Exp $', sets, pl);
392 ii.setModifier(false);
393 ii.setOutmin(2);
394 ii.setOutmax(3);
395 end
396
397 %--------------------------------------------------------------------------
398 % Get Default Plist
399 %--------------------------------------------------------------------------
400
401 function plout = getDefaultPlist()
402 persistent pl;
403 if exist('pl', 'var')==0 || isempty(pl)
404 pl = buildplist();
405 end
406 plout = pl;
407 end
408
409
410 function pl = buildplist()
411
412
413 pl = ao.getInfo('getdof', 'Default').plists;
414
415 % Conf
416 p = param({'Conf', ['Required percentage confidence level.<br>' ...
417 'It is a number between 0 and 100.']}, ...
418 {1, {95}, paramValue.OPTIONAL});
419 pl.pset(p);
420
421 % DOF
422 p = param({'dof', ['Degrees of freedom of the estimator. If it is<br>'...
423 'left empty they are calculated.']}, paramValue.EMPTY_DOUBLE);
424 pl.pset(p);
425 end
426 % END
427