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.array([[1,2,3],[2,1,3],[1,1,1]], dtype='d')
>>> 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

der(x)[source]

Evaluate the polynomial derivative at x.

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

plot(ci=1, pi=0)[source]
rnorm = None

2-norm of the residuals $$||y-Ax||_2$$

std

solution standard deviation

var

solution variance