Hi,
the attached patch does some corrections and code-optimization to
optim/inst/leasqr.m.
Except minor changes to the comment blocks, all changes apply to the
code after the actual parameter estimation.
Comments on changes which alter the output:
To get the covariance matrices of the residuals and of the data, a
certain matrix must be divided by the number of residuals "m" or the
number of degrees of freedom "m - n", respectively. The previous code
erroneously divided by "m - n" for the cov.-matr. of residuals "covr",
and nevertheless calculated the cov.-matr. of data "Vy" from the
latter by multiplying with "m / (m - n)" (reduced to "1/(1-n/m)")), so
"Vy" was also not correct. The two covariance matrices are now
calculated correctly, although of the latter only the inverse is
computed due to some code-optimization. Affected output: "covr",
"covp", "corp".
Although I found almost no reference for this, I think it is more
sensible to normalize the standard residuals "stdresid" against the
square root of the variance of the _residuals_, not the sqare root of
the variance of the _data_. This has been changed accordingly. The
former normalization was also wrong because it used the wrongly
calculated "Vy" (see above). Affected output: "stdresid"
all changes:
- corrected calculation of covariance matrix of residuals and
covariance matrix of data, thereby correcting covariance matrix of
parameters and correlation matrix of parameters
- simplified calculation of covariance matrix of parameters
- standard residuals now normalized against square root of variance of
residuals, not of data
- some repetitive computation of the same data avoided in calculation
of covariance matrix of residuals and later, one repetitive
multiplication of a matrix with a scalar avoided
- removed abs() on non-negative data in line 376
- changed hint to article of Delforge to apply to the currently
existing variables
- helptext: added comment on bounds
Except code-optimization, the above numeric changes can be
equivalently achieved simply with:
--- leasqr.m 2009-12-16 13:13:57.000000000 +0100
+++ tp_leasqr.m 2009-12-16 13:23:09.000000000 +0100
@@ -375,7 +375,7 @@
Qinv=diag(wt.*wt);
end
resid=y-f; %un-weighted residuals
- covr=resid'*Qinv*resid*Q/(m-n); %covariance of residuals
+ covr=resid'*Qinv*resid*Q/m; %covariance of residuals
Vy=1/(1-n/m)*covr; % Eq. 7-13-22, Bard %covariance of the data
jtgjinv=inv(jac'*Qinv*jac); %argument of inv may be
@@ -387,11 +387,11 @@
if (exist('sparse'))
covr=spdiags(covr,0);
- stdresid=resid./sqrt(spdiags(Vy,0));
+ stdresid=resid./sqrt(covr);
else
covr=diag(covr); % convert returned values to
% compact storage
- stdresid=resid./sqrt(diag(Vy)); % compute then convert for compact storage
+ stdresid=resid./sqrt(covr); % compute then convert for compact storage
end
Z=((m-n)*jac'*Qinv*jac)/(n*resid'*Qinv*resid);
You could compare this with the full patch to test whether the same
results are returned (they are, up to 12 digits for "covp" in my
test).
Regards, Olaf
Index: leasqr.m
===================================================================
--- leasqr.m (revision 6652)
+++ leasqr.m (working copy)
@@ -74,6 +74,9 @@
%% for each parameter. Default: [-Inf, Inf] in each row. If this
%% field is used with an existing user-side function for 'dFdp'
%% (see above) the functions interface might have to be changed.
+ %% _Warning_: If bounds are set, returned guesses of corp, covp, and
+ %% Z are generally invalid, even if final parameters are not at the
+ %% bounds.
%%
%% OUTPUT VARIABLES
%% f = column vector of values computed: f = F(x,p).
@@ -89,6 +92,10 @@
%%
%% All Zero guesses not acceptable
+ %% The following two blocks of comments are chiefly from the original
+ %% version for Matlab. For later changes the logs of the Octave Forge
+ %% svn repository should also be consulted.
+
%% A modified version of Levenberg-Marquardt
%% Non-Linear Regression program previously submitted by R.Schrager.
%% This version corrects an error in that version and also provides
@@ -120,6 +127,7 @@
%% 5) revise estimates of corp, stdev
%% Modified by Ray Muzic 11-Oct-1992
%% 1) revise estimate of Vy. remove chi2, add Z as return values
+ %% (later remark: the current code contains no variable Vy)
%% Modified by Ray Muzic 7-Jan-1994
%% 1) Replace ones(x) with a construct that is compatible with versions
%% newer and older than v 4.1.
@@ -367,37 +375,42 @@
%% diag(1/wt.^2).
%% cov matrix of data est. from Bard Eq. 7-5-13, and Row 1 Table 5.1
+ tp = wt.^2;
if (exist('sparse')) % save memory
- Q=sparse(1:m,1:m,1./wt.^2);
- Qinv=sparse(1:m,1:m,wt.^2);
+ Q = sparse (1:m, 1:m, 1 ./ tp);
+ Qinv = sparse (1:m, 1:m, tp);
else
- Q=diag((0*wt+1)./(wt.^2));
- Qinv=diag(wt.*wt);
+ Q = diag (ones (m, 1) ./ tp);
+ Qinv = diag (tp);
end
resid=y-f; %un-weighted residuals
- covr=resid'*Qinv*resid*Q/(m-n); %covariance of residuals
- Vy=1/(1-n/m)*covr; % Eq. 7-13-22, Bard %covariance of the data
+ tp = resid.' * Qinv * resid;
+ covr = (tp / m) * Q; %covariance of residuals
- jtgjinv=inv(jac'*Qinv*jac); %argument of inv may be
- %singular
- covp=jtgjinv*jac'*Qinv*Vy*Qinv*jac*jtgjinv; % Eq. 7-5-13, Bard %cov of
+ Qinv = ((m - n) / tp) * Qinv; % Qinv now contains the inverse of the
+ % guessed covariance matrix of the data.
+ % No new variable was used, but later
+ % calculations and comments using Qinv
+ % remain valid with the new Qinv.
+ %% argument of inv may be singular
+ covp = inv (jac.' * Qinv * jac); % simplified Eq. 7-5-13, Bard %cov of
% parm est
- d=sqrt(abs(diag(covp)));
+ d=sqrt(diag(covp));
corp=covp./(d*d');
if (exist('sparse'))
covr=spdiags(covr,0);
- stdresid=resid./sqrt(spdiags(Vy,0));
+ stdresid=resid ./ sqrt (covr);
else
covr=diag(covr); % convert returned values to
% compact storage
- stdresid=resid./sqrt(diag(Vy)); % compute then convert for compact storage
+ stdresid=resid ./ sqrt (covr);
end
Z=((m-n)*jac'*Qinv*jac)/(n*resid'*Qinv*resid);
%%% alt. est. of cov. mat. of parm.:(Delforge, Circulation, 82:1494-1504, 1990
%%disp('Alternate estimate of cov. of param. est.')
- %%acovp=resid'*Qinv*resid/(m-n)*jtgjinv
+ %%acovp=resid'*Qinv*resid/(m-n)*inv(jac'*Qinv*jac);
%%Calculate R^2 (Ref Draper & Smith p.46)
%%
------------------------------------------------------------------------------
This SF.Net email is sponsored by the Verizon Developer Community
Take advantage of Verizon's best-in-class app development support
A streamlined, 14 day to market process makes app distribution fast and easy
Join now and get one step closer to millions of Verizon customers
http://p.sf.net/sfu/verizon-dev2dev
_______________________________________________
Octave-dev mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/octave-dev
