octave_packages/optim-1.2.0/wsolve.m

   1 ## Author: Paul Kienzle <pkienzle@gmail.com>
   2 ## This program is granted to the public domain.
   3
   4 ## [x,s] = wsolve(A,y,dy)
   5 ##
   6 ## Solve a potentially over-determined system with uncertainty in
   7 ## the values.
   8 ##
   9 ##     A x = y +/- dy
  10 ##
  11 ## Use QR decomposition for increased accuracy.  Estimate the
  12 ## uncertainty for the solution from the scatter in the data.
  13 ##
  14 ## The returned structure s contains
  15 ##
  16 ##    normr = sqrt( A x - y ), weighted by dy
  17 ##    R such that R'R = A'A
  18 ##    df = n-p, n = rows of A, p = columns of A
  19 ##
  20 ## See polyconf for details on how to use s to compute dy.
  21 ## The covariance matrix is inv(R'*R).  If you know that the
  22 ## parameters are independent, then uncertainty is given by
  23 ## the diagonal of the covariance matrix, or
  24 ##
  25 ##    dx = sqrt(N*sumsq(inv(s.R'))')
  26 ##
  27 ## where N = normr^2/df, or N = 1 if df = 0.
  28 ##
  29 ## Example 1: weighted system
  30 ##
  31 ##    A=[1,2,3;2,1,3;1,1,1]; xin=[1;2;3];
  32 ##    dy=[0.2;0.01;0.1]; y=A*xin+randn(size(dy)).*dy;
  33 ##    [x,s] = wsolve(A,y,dy);
  34 ##    dx = sqrt(sumsq(inv(s.R'))');
  35 ##    res = [xin, x, dx]
  36 ##
  37 ## Example 2: weighted overdetermined system  y = x1 + 2*x2 + 3*x3 + e
  38 ##
  39 ##    A = fullfact([3,3,3]); xin=[1;2;3];
  40 ##    y = A*xin; dy = rand(size(y))/50; y+=dy.*randn(size(y));
  41 ##    [x,s] = wsolve(A,y,dy);
  42 ##    dx = s.normr*sqrt(sumsq(inv(s.R'))'/s.df);
  43 ##    res = [xin, x, dx]
  44 ##
  45 ## Note there is a counter-intuitive result that scaling the
  46 ## uncertainty in the data does not affect the uncertainty in
  47 ## the fit.  Indeed, if you perform a monte carlo simulation
  48 ## with x,y datasets selected from a normal distribution centered
  49 ## on y with width 10*dy instead of dy you will see that the
  50 ## variance in the parameters indeed increases by a factor of 100.
  51 ## However, if the error bars really do increase by a factor of 10
  52 ## you should expect a corresponding increase in the scatter of
  53 ## the data, which will increase the variance computed by the fit.
  54
  55 function [x_out,s]=wsolve(A,y,dy)
  56   if nargin < 2, usage("[x dx] = wsolve(A,y[,dy])"); end
  57   if nargin < 3, dy = []; end
  58
  59   [nr,nc] = size(A);
  60   if nc > nr, error("underdetermined system"); end
  61
  62   ## apply weighting term, if it was given
  63   if prod(size(dy))==1
  64     A = A ./ dy;
  65     y = y ./ dy;
  66   elseif ~isempty(dy)
  67     A = A ./ (dy * ones (1, columns(A)));
  68     y = y ./ dy;
  69   endif
  70
  71   ## system solution: A x = y => x = inv(A) y
  72   ## QR decomposition has good numerical properties:
  73   ##   AP = QR, with P'P = Q'Q = I, and R upper triangular
  74   ## so
  75   ##   inv(A) y = P inv(R) inv(Q) y = P inv(R) Q' y = P (R \ (Q' y))
  76   ## Note that b is usually a vector and Q is matrix, so it will
  77   ## be faster to compute (y' Q)' than (Q' y).
  78   [Q,R,p] = qr(A,0);
  79   x = R\(y'*Q)';
  80   x(p) = x;
  81
  82   s.R = R;
  83   s.R(:,p) = R;
  84   s.df = nr-nc;
  85   s.normr = norm(y - A*x);
  86
  87   if nargout == 0,
  88     cov = s.R'*s.R
  89     if s.df, normalized_chisq = s.normr^2/s.df, end
  90     x = x'
  91   else
  92     x_out = x;
  93   endif
  94
  95 ## We can show that uncertainty dx = sumsq(inv(R'))' = sqrt(diag(inv(A'A))).
  96 ##
  97 ## Rather than calculate inv(A'A) directly, we are going to use the QR
  98 ## decomposition we have already computed:
  99 ##
 100 ##    AP = QR, with P'P = Q'Q = I, and R upper triangular
 101 ##
 102 ## so
 103 ##
 104 ##    A'A = PR'Q'QRP' = PR'RP'
 105 ##
 106 ## and
 107 ##
 108 ##    inv(A'A) = inv(PR'RP') = inv(P')inv(R'R)inv(P) = P inv(R'R) P'
 109 ##
 110 ## For a permutation matrix P,
 111 ##
 112 ##    diag(PXP') = P diag(X)
 113 ##
 114 ## so
 115 ##    diag(inv(A'A)) = diag(P inv(R'R) P') = P diag(inv(R'R))
 116 ##
 117 ## For R upper triangular, inv(R') = inv(R)' so inv(R'R) = inv(R)inv(R)'.
 118 ## Conveniently, for X upper triangular, diag(XX') = sumsq(X')', so
 119 ##
 120 ##    diag(inv(A'A)) = P sumsq(inv(R)')'
 121 ##
 122 ## This is both faster and more accurate than computing inv(A'A)
 123 ## directly.
 124 ##
 125 ## One small problem:  if R is not square then inv(R) does not exist.
 126 ## This happens when the system is underdetermined, but in that case
 127 ## you shouldn't be using wsolve.
 128