1 ## Copyright (C) 2011 Olaf Till <olaf.till@uni-jena.de>
3 ## This program is free software; you can redistribute it and/or modify
4 ## it under the terms of the GNU General Public License as published by
5 ## the Free Software Foundation; either version 3 of the License, or
6 ## (at your option) any later version.
8 ## This program is distributed in the hope that it will be useful,
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 ## GNU General Public License for more details.
13 ## You should have received a copy of the GNU General Public License
14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
16 ## This is based on Bard, Nonlinear Parameter Estimation, Academic
17 ## Press, 1974, section 7-5, but also on own re-calculations for the
18 ## specific case of weighted least squares. The part with only certain
19 ## elements of covp or corp being defined is self-made, I don't know a
22 function hook = __covp_corp_wls__ (hook)
24 ## compute jacobian, if not already done
25 if (! isfield (hook, "jac"))
26 hook.jac = hook.dfdp (hook.pfin, hook);
28 jact = (jac = hook.jac).';
30 ## compute guessed covariance matrix of data, if not already done
31 if (! isfield (hook, "covd"))
32 hook = hook.funs.covd (hook);
34 covd_inv = inv (covd = hook.covd);
36 if (rcond (A = jact * covd_inv * jac) > eps)
38 covp = hook.covp = inv (A);
40 d = sqrt (diag (covp));
42 hook.corp = covp ./ (d * d.');
50 ## Now we have the equation "A * covp * A.' == A".
52 ## Find a particular solution for "covp * A.'".
55 ## Find a particular solution for "covp". Only uniquely defined
56 ## elements (identified later) will be further used.
57 part_covp = A \ part_covp_At.';
59 ## Find a basis for the nullspace of A.
60 if (true) # test for Octave version once submitted patch is applied
61 # to Octave (bug #33503)
62 null = @ __null_optim__;
64 if (isempty (basis = null (A)))
65 error ("internal error, singularity assumed, but null-space computed to be zero-dimensional");
68 ## Find an index (applied to both row and column) of uniquely
69 ## defined elements of covp.
70 idun = all (basis == 0, 2);
72 ## Fill in these elements.
73 covp(idun, idun) = part_covp(idun, idun);
75 ## Compute corp as far as possible at the moment.
76 d = sqrt (diag (covp));
77 corp = covp ./ (d * d.');
79 ## All diagonal elements of corp should be one, even those as yet
81 corp(1 : n + 1 : n * n) = 1;
83 ## If there are indices, applied to both row and column, so that
84 ## indexed elements within one row or one column are determined up
85 ## to a multiple of a vector, find both these vectors and the
86 ## respective indices. In the same run, use them to further fill in
87 ## corp as described below.
88 for id = 1 : (cb = columns (basis))
90 all (basis(:, (1 : cb) != id) == 0, 2) & \
92 vec = sign (basis(idx, id));
94 ## Depending on "vec", single coefficients of correlation
95 ## indexed by "idx" are either +1 or -1.
96 corp(idx, idx) = vec * vec.';