[SOLVED] Python: uncertainties and scipy.optimize.leastsq

[Link]

I've recently found out about the uncertainties package for Python, and I have been using it together with SciPy to crunch some data. Initially, the complexity of my analysis was limited to solving a simple scalar-valued scalar function, so I could just wrap fsolve. However, I'm now trying to do something more complex, and I can't figure out how to do this in a way that makes it possible to use the uncertainties package.

The situation is as follows: I have a system of equations f(...)=g(...) which are not necessarily analytically solvable. In particular, I have three equations, with a total of three known parameters and three unknowns. Let x1, x2 and x3 be my unknowns and a1, a2 and a3 my knowns, then my current code is something like this:

Code: Select all

def fitfunc(x, a1, a2, a3):
    x1, x2, x3 = x
    # ...
    return (f1(x1, x2, x3, a1, a2, a3)-g1(x1, x2, x3, a1, a2, a3),
        f2(x1, x2, x3, a1, a2, a3)-g2(x1, x2, x3, a1, a2, a3),
        f2(x1, x2, x3, a1, a2, a3)-g3(x1, x2, x3, a1, a2, a3))

x1, x2, x3 = scipy.optimize.leastsq(fitfunc, (x1_guess, x2_guess, x3_guess), args=(a1, a2, a3))[0]

What I'd like to do is use the uncertainties package to attach uncertainties to a1, a2 and a3 (i.e. use the uncertainties.Variable class for them), and get back x1, x2 and x3 with propagated uncertainties. Unfortunately, as far as I can tell, I need a function that returns a single float in order to use the wrap function from the uncertainties package. Is there any way to do what I want, or do I have to use the Monte Carlo method for this?

ETA: figured it out. You have to create a wrapper for each index of the returned tuple. Something like

Code: Select all

gen_fitfunc = lambda idx: lambda a1, a2, a3, x1_guess, x2_guess, x3_guess: \
        scipy.optimize.leastsq(fitfunc, (x1_guess, x2_guess, x3_guess), args=(a1, a2, a3))[0][idx]
x1 = uncrt.wrap(gen_fitfunc(0))(a1, a2, a3, x1_guess, x2_guess, x3_guess)
x2 = uncrt.wrap(gen_fitfunc(1))(a1, a2, a3, x1_guess, x2_guess, x3_guess)
x3 = uncrt.wrap(gen_fitfunc(2))(a1, a2, a3, x1_guess, x2_guess, x3_guess)

This means you have to call leastsq once for every parameter you want to solve for, but it should be possible to set up a cache if things get too slow.