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
|
<html lang="en">
<head>
<title>wsolve - optim_doc</title>
<meta http-equiv="Content-Type" content="text/html">
<meta name="description" content="optim_doc">
<meta name="generator" content="makeinfo 4.13">
<link title="Top" rel="start" href="index.html#Top">
<link rel="up" href="Residual-optimization.html#Residual-optimization" title="Residual optimization">
<link rel="prev" href="LinearRegression.html#LinearRegression" title="LinearRegression">
<link rel="next" href="linprog.html#linprog" title="linprog">
<link href="http://www.gnu.org/software/texinfo/" rel="generator-home" title="Texinfo Homepage">
<!--
Additional documentation for the optim package for Octave.
Copyright (C) <Olaf Till <i7tiol@t-online.de>>
You can redistribute this documentation and/or modify it under the terms
of the GNU General Public License as published by the Free Software
Foundation; either version 3 of the License, or (at your option) any
later version.
This documentation is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
Public License for more details.
You should have received a copy of the GNU General Public License along
with this documentation; if not, see <http://www.gnu.org/licenses/>.-->
<meta http-equiv="Content-Style-Type" content="text/css">
<style type="text/css"><!--
pre.display { font-family:inherit }
pre.format { font-family:inherit }
pre.smalldisplay { font-family:inherit; font-size:smaller }
pre.smallformat { font-family:inherit; font-size:smaller }
pre.smallexample { font-size:smaller }
pre.smalllisp { font-size:smaller }
span.sc { font-variant:small-caps }
span.roman { font-family:serif; font-weight:normal; }
span.sansserif { font-family:sans-serif; font-weight:normal; }
--></style>
</head>
<body>
<div class="node">
<a name="wsolve"></a>
<p>
Next: <a rel="next" accesskey="n" href="linprog.html#linprog">linprog</a>,
Previous: <a rel="previous" accesskey="p" href="LinearRegression.html#LinearRegression">LinearRegression</a>,
Up: <a rel="up" accesskey="u" href="Residual-optimization.html#Residual-optimization">Residual optimization</a>
<hr>
</div>
<h3 class="section">2.13 Function wsolve, another linear solver</h3>
<p><a name="index-wsolve-63"></a>
<h4 class="subheading">Helptext:</h4>
<p><a name="XREFwsolve"></a>
<pre class="verbatim">[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.
</pre>
<!-- -->
</body></html>
|
