Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/svd_fit.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 % SVD_FIT estimates parameters for a linear model using SVD
2 %
3 % DESCRIPTION: SVD_FIT estimates parameters for a linear model using SVD
4 %
5 % CALL: X = svd_fit([C1 C2 ... CN], Y, pl)
6 % X = svd_fit(C1,C2,C3,...,CN, Y, pl)
7 %
8 % INPUTS: C1...CN - AOs defing the models to fit the measurement set to.
9 % Y - AO which represents the measurement set
10 %
11 % Note: the length of the vectors in Ci and Y must be the same.
12 % Note: the last input AO is taken as Y.
13 %
14 % pl - parameter list (see below)
15 %
16 % OUTPUTs: X - An AO with N elements with the fitting coefficients to y_i
17 % OR
18 % - a vector of N AOs each with one fitting coefficient to y_i
19 %
20 % The procinfo field of the output AOs is filled with the following key/value
21 % pairs:
22 %
23 % 'STDX' - standard deviations of the parameters
24 % 'MSE' - the mean-squared errors
25 % 'COV' - the covariance matrix
26 %
27 %
28 % PARAMETERS:
29 %
30 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'svd_fit')">Parameters Description</a>
31 %
32 % VERSION: $Id: svd_fit.m,v 1.6 2011/04/08 08:56:12 hewitson Exp $
33 %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 function varargout = svd_fit(varargin)
36
37 % Check if this is a call for parameters
38 if utils.helper.isinfocall(varargin{:})
39 varargout{1} = getInfo(varargin{3});
40 return
41 end
42
43 import utils.const.*
44 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
45
46 % Collect input variable names
47 in_names = cell(size(varargin));
48 for ii = 1:nargin,in_names{ii} = inputname(ii);end
49
50 % Collect all AOs and plists
51 [A, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
52 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
53
54 if nargout == 0
55 error('### svd_fit can not be used as a modifier method. Please give one output');
56 end
57
58 % combine plists
59 pl = parse(pl, getDefaultPlist());
60
61 % Build matrices for fit
62
63 C = A(1:end-1);
64 Y = A(end);
65 H = C.y;
66 y = Y.y;
67 [u,s,v] = svd(H,0);
68 P = v/s*u'*y;
69 f = zeros(length(H),1); %y = zeros(length(d),1);
70 for kk = 1:length(P)
71 f = f + P(kk).*H(:,kk);
72 end
73 MSE = sum(abs(y-f).^2)./length(y);
74 a = H'*H;
75 S = inv(a)*MSE;
76 STDX = sqrt(diag(S));
77
78 % Build X
79 if find(pl,'vector_out')
80 for jj = 1:length(P)
81 X(jj) = ao(P(jj));
82 X(jj).data.setYunits(Y.yunits/C(jj).yunits);
83 X(jj).data.setDy(STDX(jj));
84 X(jj).name = sprintf('svd_fit(%s)', Y.name);
85 X(jj).addHistory(getInfo('None'), pl, ao_invars, [A(:).hist]);
86 % Set proc info
87 X(jj).procinfo = plist('STDX', STDX(jj), 'MSE', MSE, 'COV', S);
88 end
89 else
90 X = ao(P);
91 X.data.setYunits(Y.yunits/C(1).yunits);
92 X.data.setDy(STDX);
93 X.name = sprintf('svd_fit(%s)', Y.name);
94 X.addHistory(getInfo('None'), pl, ao_invars, [A(:).hist]);
95 % Set proc info
96 X.procinfo = plist('STDX', STDX, 'MSE', MSE, 'COV', S);
97 end
98
99 % Set outputs
100 varargout{1} = X;
101
102 end
103
104 %--------------------------------------------------------------------------
105 % Get Info Object
106 %--------------------------------------------------------------------------
107 function ii = getInfo(varargin)
108 if nargin == 1 && strcmpi(varargin{1}, 'None')
109 sets = {};
110 pl = [];
111 else
112 sets = {'Default'};
113 pl = getDefaultPlist;
114 end
115 % Build info object
116 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.op, '$Id: svd_fit.m,v 1.6 2011/04/08 08:56:12 hewitson Exp $', sets, pl);
117 end
118
119 %--------------------------------------------------------------------------
120 % Get Default Plist
121 %--------------------------------------------------------------------------
122 function plout = getDefaultPlist()
123 persistent pl;
124 if exist('pl', 'var')==0 || isempty(pl)
125 pl = buildplist();
126 end
127 plout = pl;
128 end
129
130 function pl = buildplist()
131 pl = plist();
132
133 % Vector out
134 p = param({'vector_out','The estimated coefficients are output as a vector of AOs.'}, paramValue.TRUE_FALSE);
135 pl.append(p);
136
137 end
138 % END