wsolve - Weighted linear and polynomial solver with uncertainty¶
wsolve |
Given a linear system \(y = A x + \delta y\), estimates \(x\) and \(\delta x\). |
wpolyfit |
Return the polynomial of degree \(n\) that minimizes \(\sum(p(x_i) - y_i)^2/\sigma_i^2\). |
LinearModel |
Model evaluator for linear solution to \(Ax = y\). |
PolynomialModel |
Model evaluator for best fit polynomial \(p(x) = y +/- \delta y\). |
Weighted linear and polynomial solver with uncertainty.
Given \(A \bar x = \bar y \pm \delta \bar y\), solve using s = wsolve(A,y,dy)
wsolve uses the singular value decomposition for increased accuracy.
The uncertainty in the solution is estimated from the scatter in the data. Estimates the uncertainty for the solution from the scatter in the data.
The returned model object s provides:
s.x solution s.std uncertainty estimate assuming no correlation s.rnorm residual norm s.DoF degrees of freedom s.cov covariance matrix s.ci(p) confidence intervals at point p s.pi(p) prediction intervals at point p s(p) predicted value at point p
Example¶
Weighted system:
>>> import numpy as np
>>> from bumps import wsolve
>>> A = np.matrix("1,2,3;2,1,3;1,1,1",'d').A
>>> dy = [0.2,0.01,0.1]
>>> y = [ 14.16, 13.01, 6.15]
>>> s = wsolve.wsolve(A,y,dy)
>>> print(", ".join("%0.2f +/- %0.2f"%(a,b) for a,b in zip(s.x,s.std)))
1.05 +/- 0.17, 2.20 +/- 0.12, 2.91 +/- 0.12
Note there is a counter-intuitive result that scaling the estimated uncertainty in the data does not affect the computed uncertainty in the fit. This is the correct result — if the data were indeed selected from a process with ten times the uncertainty, you would expect the scatter in the data to increase by a factor of ten as well. When this new data set is fitted, it will show a computed uncertainty increased by the same factor. Monte carlo simulations bear this out. The conclusion is that the dataset carries its own information about the variance in the data, and the weight vector serves only to provide relative weighting between the points.
-
bumps.wsolve.
wsolve
(A, y, dy=1, rcond=1e-12)[source]¶ Given a linear system \(y = A x + \delta y\), estimates \(x\) and \(\delta x\).
A is an n x m array of measurement points.
y is an n x k array or vector of length n of measured values at A.
dy is a scalar or an n x 1 array of uncertainties in the values at A.
Returns
LinearModel
.
-
bumps.wsolve.
wpolyfit
(x, y, dy=1, degree=None, origin=False)[source]¶ Return the polynomial of degree \(n\) that minimizes \(\sum(p(x_i) - y_i)^2/\sigma_i^2\).
if origin is True, the fit should go through the origin.
Returns
PolynomialModel
.
-
class
bumps.wsolve.
LinearModel
(x=None, DoF=None, SVinv=None, rnorm=None)[source]¶ Bases:
object
Model evaluator for linear solution to \(Ax = y\).
Use s(A) to compute the predicted value of the linear model s at points given on the rows of \(A\).
Computes a confidence interval (range of likely values for the mean at \(x\)) or a prediction interval (range of likely values seen when measuring at \(x\)). The prediction interval gives the width of the distribution at \(x\). This should be the same regardless of the number of measurements you have for the value at \(x\). The confidence interval gives the uncertainty in the mean at \(x\). It should get smaller as you increase the number of measurements. Error bars in the physical sciences usually show a \(1-\alpha\) confidence value of \(\text{erfc}(1/\sqrt{2})\), representing a \(1-\sigma\) standand deviation of uncertainty in the mean.
Confidence intervals for the expected value of the linear system evaluated at a new point \(w\) are given by the \(t\) distribution for the selected interval \(1-\alpha\), the solution \(x\), and the number of degrees of freedom \(n-p\):
\[w^T x \pm t^{\alpha/2}_{n-p} \sqrt{ \text{var}(w) }\]where the variance \(\text{var}(w)\) is given by:
\[\text{var}(w) = \sigma^2 (w^T (A^TA)^{-1} w)\]Prediction intervals are similar, except the variance term increases to include both the uncertainty in the predicted value and the variance in the data:
\[\text{var}(w) = \sigma^2 (1 + w^T (A^TA)^{-1} w)\]-
DoF
= None¶ number of degrees of freedom in the solution
-
ci
(A, sigma=1)[source]¶ Compute the calculated values and the confidence intervals for the linear model evaluated at \(A\).
sigma=1 corresponds to a \(1-\sigma\) confidence interval
Confidence intervals are sometimes expressed as \(1-\alpha\) values, where \(\alpha = \text{erfc}(\sigma/\sqrt{2})\).
-
cov
¶ covariance matrix [inv(A’A); O(n^3)]
-
p
¶ p-value probability of rejection
-
pi
(A, p=0.05)[source]¶ Compute the calculated values and the prediction intervals for the linear model evaluated at \(A\).
p=0.05 corresponds to the 95% prediction interval.
-
rnorm
= None¶ 2-norm of the residuals \(||y-Ax||_2\)
-
std
¶ solution standard deviation [sqrt(var); O(n^2)]
-
var
¶ solution variance [diag(cov); O(n^2)]
-
x
= None¶ solution to the equation \(Ax = y\)
-
-
class
bumps.wsolve.
PolynomialModel
(x, y, dy, s, origin=False)[source]¶ Bases:
object
Model evaluator for best fit polynomial \(p(x) = y +/- \delta y\).
Use p(x) for PolynomialModel p to evaluate the polynomial at all points in the vector x.
-
DoF
= None¶ number of degrees of freedom in the solution
-
ci
(x, sigma=1)[source]¶ Evaluate the polynomial and the confidence intervals at x.
sigma=1 corresponds to a 1-sigma confidence interval
-
coeff
= None¶ polynomial coefficients
-
cov
¶ covariance matrix
Note that the ones column will be absent if origin is True.
-
degree
= None¶ polynomial degree
-
origin
= None¶ True if polynomial goes through the origin
-
p
¶ p-value probability of rejection
-
pi
(x, p=0.05)[source]¶ Evaluate the polynomial and the prediction intervals at x.
p = 1-alpha = 0.05 corresponds to 95% prediction interval
-
rnorm
= None¶ 2-norm of the residuals \(||y-Ax||_2\)
-
std
¶ solution standard deviation
-
var
¶ solution variance
-