1 ## Author: Paul Kienzle <pkienzle@gmail.com>
2 ## This program is granted to the public domain.
4 ## [x,s] = wsolve(A,y,dy)
6 ## Solve a potentially over-determined system with uncertainty in
11 ## Use QR decomposition for increased accuracy. Estimate the
12 ## uncertainty for the solution from the scatter in the data.
14 ## The returned structure s contains
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
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
25 ## dx = sqrt(N*sumsq(inv(s.R'))')
27 ## where N = normr^2/df, or N = 1 if df = 0.
29 ## Example 1: weighted system
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'))');
37 ## Example 2: weighted overdetermined system y = x1 + 2*x2 + 3*x3 + e
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);
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.
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
60 if nc > nr, error("underdetermined system"); end
62 ## apply weighting term, if it was given
67 A = A ./ (dy * ones (1, columns(A)));
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
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).
85 s.normr = norm(y - A*x);
89 if s.df, normalized_chisq = s.normr^2/s.df, end
95 ## We can show that uncertainty dx = sumsq(inv(R'))' = sqrt(diag(inv(A'A))).
97 ## Rather than calculate inv(A'A) directly, we are going to use the QR
98 ## decomposition we have already computed:
100 ## AP = QR, with P'P = Q'Q = I, and R upper triangular
104 ## A'A = PR'Q'QRP' = PR'RP'
108 ## inv(A'A) = inv(PR'RP') = inv(P')inv(R'R)inv(P) = P inv(R'R) P'
110 ## For a permutation matrix P,
112 ## diag(PXP') = P diag(X)
115 ## diag(inv(A'A)) = diag(P inv(R'R) P') = P diag(inv(R'R))
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
120 ## diag(inv(A'A)) = P sumsq(inv(R)')'
122 ## This is both faster and more accurate than computing inv(A'A)
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.