% [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.