1

Given a target function of this type:

<a href="http://www.codecogs.com/eqnedit.php?latex=\dpi{200}&space;y^{(i)}&space;=&space;\sum_{j=1}^{N_1}&space;\sum_{k=1}^{N_2}&space;D_j*\left&space;[&space;[\exp(-a_{j}*x_{jk}^{(i)})-1&space;]^2&space;-&space;1&space;\right&space;]" target="_blank"><img src="http://latex.codecogs.com/gif.latex?\dpi{200}&space;y^{(i)}&space;=&space;\sum_{j=1}^{N_1}&space;\sum_{k=1}^{N_2}&space;D_j*\left&space;[&space;[\exp(-a_{j}*x_{jk}^{(i)})-1&space;]^2&space;-&space;1&space;\right&space;]" title="y^{(i)} = \sum_{j=1}^{N_1} \sum_{k=1}^{N_2} D_j*\left [ [\exp(-a_{j}*x_{jk}^{(i)})-1 ]^2 - 1 \right ]"
  /></a>

where D_j and a_j are the parameters and the summations over index j, k are not fixed in number (N_1 and N_2 may vary). For a set of input data (x, y) (x is a 2D matrix), how to fit the parameters involved using numpy, scipy, or lmfit-py (https://github.com/lmfit/lmfit-py/blob/master/doc/intro.rst)?

There is a relevant post here Fitting a sum to data in Python, but my case seems to be a bit more complicated. Thanks for any comment!

Huang Bing
  • 13
  • 4
  • 4
    can you make a readable formula elsewhere, e.g. word/excel and put a screenshot of the formula here? would be better. Thanks. – Rockbar Nov 17 '17 at 20:36
  • If I understand correctly, your summation range is arbitrary, hence you have an infinite amount of D_j and a_j to optimize the fit. This is not a well-defined problem with a definite solution. – Christoph Nov 17 '17 at 21:09
  • @Rockbar Thanks for the suggestion, now you can see the reasable formula by pressing the button `Run code snippet. :) – Huang Bing Nov 17 '17 at 21:20
  • @Christoph Thanks for the prompt reply! Now I updated the formula and make the form more explicit, i.e., there the number of parameters D_j and a_j are finite. – Huang Bing Nov 17 '17 at 21:22
  • You still have an infinite amount of degrees of freedom when N1 and N2 are not fixed – Christoph Nov 17 '17 at 21:39
  • @Christoph ok, say let's fix N1 to 4, N2 to 5; if in some case `i, `k runs from 1 to 3, then we can pad the 4th and 5th column of **x** by zero. Now it should be ok. – Huang Bing Nov 17 '17 at 22:51
  • Why do you think that the solution you linked to is not valid for your case, only for your double sum? – mikuszefski Nov 20 '17 at 08:42
  • @mikuszefski yes, you are right...I actually made it work after posting this messages. Nevertheless, the answer by M Newville below is more elegent! – Huang Bing Nov 21 '17 at 11:56

1 Answers1

0

I'm not sure your formula is enough to give a complete answer, but if I understand correctly, then you will know N_1 before the fit, it will be some finite number, and you really want 2*N_1 fit parameters. If that is correct, then what you want to do is generate and use 2*N_1 parameters dynamically. To do that, you can set up the objective function something like this:

def objective(params, ydata, xdata, n_1):
    npts = len(ydata)
    ymodel = np.zeros(npts)
    for j in range(n_1):
        dj = params['d_%i' % (i+1)].value
        aj = params['a_%i' % (i+1)].value
        submodel = calc_model(dj, aj, xdata)
        ymodel += submodel

    return (ymodel - ydata)  

And build the Parameters and run the fit like this:

 params = Parameters()
 for i in range(n_1):
     params.add('d_%i' % (i+1), value=1.0)  # set init values here
     params.add('a_%i' % (i+1), value=0.2)  # 

result = minimize(objective, params, args=(ydata, xdata, n_1))
M Newville
  • 7,486
  • 2
  • 16
  • 29