Fit simulated and experimental data points with Python

For fitting your data with an arbitrary function you usually want the Levenberg–Marquardt algorithm, and that is what scipy.optimize.leastsq uses, so you are most likely using the correct function.

Have a look at the tutorial in section 5.4 of this page and see if it helps.

It is also possible your underlying function is difficult to fit (what is the function?), but it sounds like you may be having other issues.

Also - as great as StackOverflow is, you'll likely get much more knowledgeable responses for scipy questions by posting directly to the Scipy-User mailing list with some sample code and more detail.

If you don't know the expected functional form of the global, but can predict an expected value for the "next" point given the current state of the system you might consider using a Kalman filter (yeah, "filter" sounds goofy in a fitting context but the name is historical and can't easily be changed now).

The underlying math can look a bit scary, but this important point is that you don't have to understand it. You typically need to be able to

1. Define a representation space
2. Express your data (simulated or experimental) in the representation space.
3. Define a update procedure that gets a "next" representation from a given one
4. Extract the desired curve from a series representation as returned by the fitter

There seems to be at least one existing python package to support this (note that the interface here is different from the ones that I'm used to and I can't offer much advice on using it).

Not exactly an answer but there is also PyMinuit to try.

http://code.google.com/p/pyminuit/

What you want to do is convert your pdf and data points to chi^2 or -ln(likelihood) or a metric of your choosing and use PyMinuit to minimize that metric. It can be configured to be very verbose so you can find out where things went wrong(if it did go wrong).

Since you only have two parameters, you should just do a grid search.

results = {}
for p0 in parameter_space_for_first_parameter:
     for p1 in parameter_space_for_second_parameter:
           results[p0,p1] = compare(p0,p1)

If you can afford the computation, compare should perform multiple runs (with different initialisations) and compute mean and standard deviations. You can try using my package jug to manage your computations (it was designed for exactly this sort of thing).

Finally, plot the results, look at the minima (there might be several). This is a "stupid" method, but it will work in many situations where other methods get stuck.

If this is too much computation, you can do this in two passes: a coarse-grained grid followed by a fine-grained one near the minima of the coarse-grained space.