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