doc: Fitting with multiple curves

[pbeaucage]

@ktoth17 was trying to do a multi-curve fit, where a power law slope was fixed between two different datasets. His example code:

def exp_func(x, para):
    A, B, C = para
    return A * x**B + C
# Error function taking an array of parameters (para), x and y
def err(para,x,y):
    return abs(exp_func(x, para)-y)
# Global error to minimize
# Parameters here are [A1 A2 B C2] since B is held constant and C1 = 0.
def err_global(para, x, y1, y2):
    p1 = para[0], para[2], 0
    p2 = para[1], para[2], para[3]
    err1 = err(p1,x,y1)
    err2 = err(p2,x,y2)
    return np.concatenate((err1,err2))
# Initial guess for parameters [A1 A2 B C2]
init_guess = [1, 1, -4, 0]
# Minimize the sum of squares of a set of equations
para_best, ier = scipy.optimize.leastsq(err_global, init_guess, args=(x_values,y1_values,y2_values))

He asked what the best way to do this was in PyHyper, a good question and an area that definitely should be in the documentation.

[pbeaucage]

Copying Slack discussion for preservation:

15 replies

Kristof Toth 4 days ago
Here I'm fitting two data sets (x,y1) and (x,y2) to a simple power law of A*x**B + C. I would like to keep the term B constant and find the optimal A, B, and C values using scipy.optimize.leastsq. Is there a way to incorporate a loop or input a tuple to have an arbitrary number of (x,y's) data sets? Thank you!

Peter Beaucage 4 days ago
So this works on a single i1 vs q / i2 vs q pair and you want to do it on an arbitrary stack of i1 vs q / i2 vs q at each of E= e1 e 2 e 3 etc?

Peter Beaucage 4 days ago
Do you want to pin the parameters between runs or float them in each energy but fix between the two i1 and i2 datasets?

Kristof Toth 4 days ago
Yes, currently it solves for the parameters [A1 A2 B C2] because B is held constant and C1 I set for 0. I would like to pin B between each intensity and q but find individual A's and C's...

Kristof Toth 4 days ago
Normally, I would think it loops or input arrays, but I haven't seen implementations where there are arbitrary error functions that get concatenated for leastsq( ). (edited)

Peter Beaucage 4 days ago
Normally the answer would be (at least for pyhyper) a split/apply/combine but I’d need to look at the xarray docs to see if you can concatenate the two arrays. One option would be to stack the xarrays along some other axis (sampleID, say) and then handle it as a single object.

Peter Beaucage 4 days ago
BUT in this case a for loop should probably work
✅
1

Peter Beaucage 4 days ago
Something like (sorry on phone):

Peter Beaucage 4 days ago
Results = {}
for en in array1.energy:
x = array1.sel(energy=en).q
y1 = array1.sel(energy=en) # may have to transform to np array here
y2 = array2.sel(energy=en) # same
Results[en] = fitfunc(x,y1,y2)

Kristof Toth 4 days ago
Ah, that's helpful!

Peter Beaucage 4 days ago
If seed parameters need to change you can pre stash them in a dict like results. And the selection queries may need to be trimmed down to reduce chi or another axis. And it assumes both scans have the same energy sets - may need to handle that edge case. Also afaik this will be lower performance than the split apply combine because it can’t automatically multithread. But it should work.

Veronica Reynolds 4 days ago
Just jumping in here to say I recently started using this lmfit package: https://lmfit.github.io/lmfit-py/index.html I've found it to be more nimble than scipy optimize for choosing what you fit or fix, setting bounds, etc. Not sure if it would have any benefit for your situation!
🙌
1

Peter Beaucage 4 days ago
Yeah, my initial instinct was that this was beyond scipy optimize but @tyler Martin corrected me. emcee and bumps are more advanced alternatives.