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