/usr/share/octave/packages/optim-1.3.0/wsolve.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 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 | ## Author: Paul Kienzle <pkienzle@gmail.com>
## This program is granted to the public domain.
## [x,s] = wsolve(A,y,dy)
##
## Solve a potentially over-determined system with uncertainty in
## the values.
##
## A x = y +/- dy
##
## Use QR decomposition for increased accuracy. Estimate the
## uncertainty for the solution from the scatter in the data.
##
## The returned structure s contains
##
## normr = sqrt( A x - y ), weighted by dy
## R such that R'R = A'A
## df = n-p, n = rows of A, p = columns of A
##
## See polyconf for details on how to use s to compute dy.
## The covariance matrix is inv(R'*R). If you know that the
## parameters are independent, then uncertainty is given by
## the diagonal of the covariance matrix, or
##
## dx = sqrt(N*sumsq(inv(s.R'))')
##
## where N = normr^2/df, or N = 1 if df = 0.
##
## Example 1: weighted system
##
## A=[1,2,3;2,1,3;1,1,1]; xin=[1;2;3];
## dy=[0.2;0.01;0.1]; y=A*xin+randn(size(dy)).*dy;
## [x,s] = wsolve(A,y,dy);
## dx = sqrt(sumsq(inv(s.R'))');
## res = [xin, x, dx]
##
## Example 2: weighted overdetermined system y = x1 + 2*x2 + 3*x3 + e
##
## A = fullfact([3,3,3]); xin=[1;2;3];
## y = A*xin; dy = rand(size(y))/50; y+=dy.*randn(size(y));
## [x,s] = wsolve(A,y,dy);
## dx = s.normr*sqrt(sumsq(inv(s.R'))'/s.df);
## res = [xin, x, dx]
##
## Note there is a counter-intuitive result that scaling the
## uncertainty in the data does not affect the uncertainty in
## the fit. Indeed, if you perform a monte carlo simulation
## with x,y datasets selected from a normal distribution centered
## on y with width 10*dy instead of dy you will see that the
## variance in the parameters indeed increases by a factor of 100.
## However, if the error bars really do increase by a factor of 10
## you should expect a corresponding increase in the scatter of
## the data, which will increase the variance computed by the fit.
function [x_out,s]=wsolve(A,y,dy)
if nargin < 2, usage("[x dx] = wsolve(A,y[,dy])"); end
if nargin < 3, dy = []; end
[nr,nc] = size(A);
if nc > nr, error("underdetermined system"); end
## apply weighting term, if it was given
if prod(size(dy))==1
A = A ./ dy;
y = y ./ dy;
elseif ~isempty(dy)
A = A ./ (dy * ones (1, columns(A)));
y = y ./ dy;
endif
## system solution: A x = y => x = inv(A) y
## QR decomposition has good numerical properties:
## AP = QR, with P'P = Q'Q = I, and R upper triangular
## so
## inv(A) y = P inv(R) inv(Q) y = P inv(R) Q' y = P (R \ (Q' y))
## Note that b is usually a vector and Q is matrix, so it will
## be faster to compute (y' Q)' than (Q' y).
[Q,R,p] = qr(A,0);
x = R\(y'*Q)';
x(p) = x;
s.R = R;
s.R(:,p) = R;
s.df = nr-nc;
s.normr = norm(y - A*x);
if nargout == 0,
cov = s.R'*s.R
if s.df, normalized_chisq = s.normr^2/s.df, end
x = x'
else
x_out = x;
endif
## We can show that uncertainty dx = sumsq(inv(R'))' = sqrt(diag(inv(A'A))).
##
## Rather than calculate inv(A'A) directly, we are going to use the QR
## decomposition we have already computed:
##
## AP = QR, with P'P = Q'Q = I, and R upper triangular
##
## so
##
## A'A = PR'Q'QRP' = PR'RP'
##
## and
##
## inv(A'A) = inv(PR'RP') = inv(P')inv(R'R)inv(P) = P inv(R'R) P'
##
## For a permutation matrix P,
##
## diag(PXP') = P diag(X)
##
## so
## diag(inv(A'A)) = diag(P inv(R'R) P') = P diag(inv(R'R))
##
## For R upper triangular, inv(R') = inv(R)' so inv(R'R) = inv(R)inv(R)'.
## Conveniently, for X upper triangular, diag(XX') = sumsq(X')', so
##
## diag(inv(A'A)) = P sumsq(inv(R)')'
##
## This is both faster and more accurate than computing inv(A'A)
## directly.
##
## One small problem: if R is not square then inv(R) does not exist.
## This happens when the system is underdetermined, but in that case
## you shouldn't be using wsolve.
|