Home > freetb4matlab > optim > wsolve.m

wsolve

PURPOSE ^

% [x,s] = wsolve(A,y,dy)

SYNOPSIS ^

function [x_out,s]=wsolve(A,y,dy)

DESCRIPTION ^

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:
Generated on Sat 16-May-2009 00:04:49 by m2html © 2003