octave_packages/optim-1.2.0/private/__covp_corp_wls__.m

   1 ## Copyright (C) 2011 Olaf Till <olaf.till@uni-jena.de>
   2 ##
   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.
   7 ##
   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.
  12 ##
  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/>.
  15
  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
  20 ## reference.
  21
  22 function hook = __covp_corp_wls__ (hook)
  23
  24   ## compute jacobian, if not already done
  25   if (! isfield (hook, "jac"))
  26     hook.jac = hook.dfdp (hook.pfin, hook);
  27   endif
  28   jact = (jac = hook.jac).';
  29
  30   ## compute guessed covariance matrix of data, if not already done
  31   if (! isfield (hook, "covd"))
  32     hook = hook.funs.covd (hook);
  33   endif
  34   covd_inv = inv (covd = hook.covd);
  35
  36   if (rcond (A = jact * covd_inv * jac) > eps)
  37
  38     covp = hook.covp = inv (A);
  39
  40     d = sqrt (diag (covp));
  41
  42     hook.corp = covp ./ (d * d.');
  43
  44   else
  45
  46     n = hook.np;
  47
  48     covp = NA (n);
  49
  50     ## Now we have the equation "A * covp * A.' == A".
  51
  52     ## Find a particular solution for "covp * A.'".
  53     part_covp_At = A \ A;
  54
  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.';
  58
  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__;
  63     endif
  64     if (isempty (basis = null (A)))
  65       error ("internal error, singularity assumed, but null-space computed to be zero-dimensional");
  66     endif
  67
  68     ## Find an index (applied to both row and column) of uniquely
  69     ## defined elements of covp.
  70     idun = all (basis == 0, 2);
  71
  72     ## Fill in these elements.
  73     covp(idun, idun) = part_covp(idun, idun);
  74
  75     ## Compute corp as far as possible at the moment.
  76     d = sqrt (diag (covp));
  77     corp = covp ./ (d * d.');
  78
  79     ## All diagonal elements of corp should be one, even those as yet
  80     ## NA.
  81     corp(1 : n + 1 : n * n) = 1;
  82
  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))
  89       if (any (idx = \
  90                all (basis(:, (1 : cb) != id) == 0, 2) & \
  91                basis(:, id) != 0))
  92         vec = sign (basis(idx, id));
  93
  94         ## Depending on "vec", single coefficients of correlation
  95         ## indexed by "idx" are either +1 or -1.
  96         corp(idx, idx) = vec * vec.';
  97       endif
  98     endfor
  99
 100     hook.covp = covp;
 101     hook.corp = corp;
 102
 103   endif
 104
 105 endfunction