solve linear equations given variables and uncertainties: scipy-optimize?

chemistry, numpy, optimization, python, scipy

Solution

This is my approach. Assuming `x1-x4` are approximately normally distributed around each mean (1-sigma uncertainty), the problem is turning into one of minimizing the sum of square of errors, with 3 linear constrain functions. Therefore, we can attack it using `scipy.optimize.fmin_slsqp()`

In [19]:

def eq_f1(x):
    return (x*np.array([0.5, 1.0, 1.5, 2.0])).sum()-8
def eq_f2(x):
    return (x*np.array([0.0, 0.0, 1.0, 1.0])).sum()-4
def eq_f3(x):
    return (x*np.array([1.0, 1.0, 0.0, 0.0])).sum()-1
def error_f(x):
    error=(x-np.array([0.246, 0.749, 1.738, 2.248]))/np.array([0.007, 0.010, 0.009, 0.007])
    return (error*error).sum()
In [20]:

so.fmin_slsqp(error_f, np.array([0.246, 0.749, 1.738, 2.248]), eqcons=[eq_f1, eq_f2, eq_f3])
Optimization terminated successfully.    (Exit mode 0)
            Current function value: 2.17576389592
            Iterations: 4
            Function evaluations: 32
            Gradient evaluations: 4
Out[20]:
array([ 0.25056582,  0.74943418,  1.74943418,  2.25056582])

Problem

I'd like to minimize a set of equations where the variables are known with their uncertainties. In essence I'd like to test the hypothesis that the given measured variables conform to the formula constraints given by the equations. This seems like something I should be able to do with scipy-optimize. For example I have three equations: ``` 8 = 0.5 * x1 + 1.0 * x2 + 1.5 * x3 + 2.0 * x4 4 = 0.0 * x1 + 0.0 * x2 + 1.0 * x3 + 1.0 * x4 1 = 1.0 * x1 + 1.0 * x2 + 0.0 * x3 + 0.0 * x4 ``` And four measured unknowns with their 1-sigma uncertainty: ``` x1 = 0.246 ± 0.007 x2 = 0.749 ± 0.010 x3 = 1.738 ± 0.009 x4 = 2.248 ± 0.007 ``` Looking for any pointers in the right direction.

Original source