How to fit an Implicit 1D function?

[Giulio de Vito]

Hi everybody,

I have this implicit function to fit in order to extract the coefficients:

 (a*y^3 + (b^2 - x^2)*y)^2 + (c*x*y)^2 -d =0

x and y are known and a,b,c,d have to be extracted by fitting the function with the y values. 

Is there a way using lmfit to extract those parameters?

It would be correct if I consider the function as multi-dimensional and calculate the residuals?

thank you

[Darshan Beniwal]

Hi Giulio de Vito 

Yes, there is a way to extract parameters using Imfit. Here I'm attaching the manual titled " Non-Linear Least-Squares Minimization and Curve-Fitting for Python Release 0.9.12". You may take help of it.

With Regards!!!      

[Matt Newville]

Hi Giulio,

Yes, this should be possible to solve with lmfit.   I think what you want to do is find a, b, c, and d such that your expression that ideally will equal zero is as small as possible in the least-squares sense (which will prevent the expression from going far below 0).

This is not a data fitting problem, so you'll want to use the `lmfit.minimize` interface.   First define an objective function to turn your parameters and extra arrays (your x and y) into the array to be minimized:

     from lmfit import Parameters, minimize, report_fit
     def obj(params, x, y):
        'objective function: returns the array to be minimized'
        a = params['a']
        b = params['b']
        c = params['c']
        d = params['d']
        return  (a*y^3 + (b^2 - x^2)*y)^2 + (c*x*y)^2 -d

Then create a Parameters object, with initial values for a, b, c, d
   
     params = Parameters()
     params.add('a', value=10.0, min=0)
     params.add('b', value=0.5) .
     params.add('c', value=0.0)
     params.add('d', value=1.0)

Then do th fit and report the result:

    result = minimizer(fcn2min, params, fcn_args=(x, data))
    report_fit(result)

Hope that's enough to get you started.

--Matt

[Giulio de Vito]

Hi everyone,

thank you for the support, I have read the guide you had linked, it was very helpful.

After some troubles I arrived to this:

def f(x, y, a, b, c, d):
return (a*y**3+(b**2-x**2)*y)**2+(y*c*x)**2-d


df = pd.read_csv(Datafile)
x= df['Frequency'][:] #Frequency along x
y= df['Amplitude'][:] #Average of amplitude


mod = lmfit.Model(f, independent_vars=['x','y']) #model declared

#Initial values and contrains for the fitting coefficients
mod.set_param_hint('a',value=1)
mod.set_param_hint('b',value=np.mean(x), min=0)
mod.set_param_hint('c',value=100, min=0)
mod.set_param_hint('d',value=1, min=0)

pars = mod.make_params(verbose=True) #Make parameters from the model

out = mod.fit(y, pars, x=x, y=y, verbose=True) #First fitting
print(mod.print_param_hints(colwidth=10)) #show the parameters
print(out.fit_report())



#Re-fitting the model

result=out.fit(y,params=out.params,x=x, y=y) #Second fitting

print(out.fit_report()) #Print results

#Evaluation of residuals

evaluation = out.eval() #Saving those numbers in name_file_eval
uncert = out.eval_uncertainty()

#calculate R-squared
lmfit_Rsquared = '%.5f'%(1 - out.residual.var()/np.var(y))

print("The R-square value is: ", lmfit_Rsquared)

print(out.ci_report())
out.params.pretty_print(colwidth=10, fmt='g')

Which gives the right results, but with the problem (I guess) that the model try to fit every single point in x and y data.

While it should just use 'x' as independent variable  and minimize the residuals with 'y'.

Thanks to Matt I also implemented the minimize function:

def obj(params, x, y):
'objective function: returns the array to be minimized'
a = params['a']
b = params['b']
c = params['c']
d = params['d']
return (a*y^3 + (b^2 - x^2)*y)^2 + (c*x*y)^2 -d

params = Parameters()
params.add('a', value=1)
params.add('b', value=np.mean(x), min=0)
params.add('c', value=100, min=0)
params.add('d', value=1, min=0)
result = minimizer(f, params, fcn_args=(x, y))
report_fit(result)

But in this case it gives the error message:

TypeError: 'module' object is not callable

There are attached the output response

[Matt Newville]

Hi Giulio,

Well, I apologize for a typo in the script I sent, but I think you also went off in sort of an unfortunate direction.

On Mon, May 13, 2019 at 5:21 AM Giulio de Vito <giulio...@gmail.com> wrote:

Hi everyone,

thank you for the support, I have read the guide you had linked, it was very helpful.

After some troubles I arrived to this:

def f(x, y, a, b, c, d):
return (a*y**3+(b**2-x**2)*y)**2+(y*c*x)**2-d


df = pd.read_csv(Datafile)
x= df['Frequency'][:] #Frequency along x
y= df['Amplitude'][:] #Average of amplitude


mod = lmfit.Model(f, independent_vars=['x','y']) #model declared

Um, but your problem is not a curve fitting problem.  Using Model for a problem that is not curve-fitting will cause many headaches, which you have seen.

#Initial values and contrains for the fitting coefficients
mod.set_param_hint('a',value=1)
mod.set_param_hint('b',value=np.mean(x), min=0)
mod.set_param_hint('c',value=100, min=0)
mod.set_param_hint('d',value=1, min=0)

Using `model.set_param_hint` is definitely anti-recommended here.   I am continually surprised to see people using this method of a model to make initial guesses for a particular fit.  It is a misunderstanding of the intent.  I think the docs are going to have to be altered to strongly discourage this.

pars = mod.make_params(verbose=True) #Make parameters from the model

out = mod.fit(y, pars, x=x, y=y, verbose=True) #First fitting
print(mod.print_param_hints(colwidth=10)) #show the parameters
print(out.fit_report())



#Re-fitting the model

result=out.fit(y,params=out.params,x=x, y=y) #Second fitting

Why are you re-doing the fit?  What is that trying to accomplish?

print(out.fit_report()) #Print results

#Evaluation of residuals

evaluation = out.eval() #Saving those numbers in name_file_eval
uncert = out.eval_uncertainty()

#calculate R-squared
lmfit_Rsquared = '%.5f'%(1 - out.residual.var()/np.var(y))

Hm, really?

print("The R-square value is: ", lmfit_Rsquared)

print(out.ci_report())
out.params.pretty_print(colwidth=10, fmt='g')

Which gives the right results, but with the problem (I guess) that the model try to fit every single point in x and y data.

While it should just use 'x' as independent variable  and minimize the residuals with 'y'.

Thanks to Matt I also implemented the minimize function:

def obj(params, x, y):
'objective function: returns the array to be minimized'
a = params['a']
b = params['b']
c = params['c']
d = params['d']
return (a*y^3 + (b^2 - x^2)*y)^2 + (c*x*y)^2 -d


params = Parameters()
params.add('a', value=1)
params.add('b', value=np.mean(x), min=0)
params.add('c', value=100, min=0)
params.add('d', value=1, min=0)
result = minimizer(f, params, fcn_args=(x, y))
report_fit(result)

But in this case it gives the error message:

TypeError: 'module' object is not callable

Well, I apologize for the typos, but it is a shame you didn't investigate that error further (or read the whole traceback) and find that it should have been 

def obj(params, x, y):

    'objective function: returns the array to be minimized'

    a = params['a']

    b = params['b']

    c = params['c']

    d = params['d']

    return (a*y**3 + (b**2 - x**2)*y)**2 + (c*x*y)**2 -d

params = Parameters()

params.add('a', value=1)

params.add('b', value=np.mean(x), min=0)

params.add('c', value=100, min=0)

params.add('d', value=1, min=0)

result = minimize(obj, params, args=(x, y))

Cheers,

--Matt 

[Giulio de Vito]

Hi Matt,

Yes as you said it could not be a problem of model fitting but rather mathematical one, to be more clear I attached an image 'Ham_plots' of what it supposes to do.

As you can see in the plot the equation  "(a*y**3 + (b**2 - x**2)*y)**2 + (c*x*y)**2 -d" gives the black points if 'a,b,c,d' are given.

But in this case the black points are given are 'a,b,c,d' have to be extracted.

One of the major challenges is that is that there are multiples "y" for one "x" so the function cannot associate uniquely the two set of data.

Moreover it is not possible to have an explicit 'y' with respect to 'x'.

In some cases the lmfit function works and it can extract parameters correctly but, in case for two consequently x (x2-x1) there is a big gap (y2-y1) The model explode and gives a 'd' value too big.

An example is in the picture 'lmfit_badresult'

Why are you re-doing the fit?  What is that trying to accomplish?

I use a second fitting to confront the first fitting with the second to make sure that the algorithm reached the steady state and there is no a significant variation respect to the initial guess.

Do you suggest to change the fitting method to avoid the explosion of the interpolating curve?

Maybe  a weighting of the data, especially near the most discontinuous region can help?

thanks,

Giulio 

[Matt Newville]

On Wed, May 15, 2019 at 4:40 AM Giulio de Vito <giulio...@gmail.com> wrote:

Hi Matt,

Yes as you said it could not be a problem of model fitting but rather mathematical one, to be more clear I attached an image 'Ham_plots' of what it supposes to do.

As you can see in the plot the equation  "(a*y**3 + (b**2 - x**2)*y)**2 + (c*x*y)**2 -d" gives the black points if 'a,b,c,d' are given.

Actually, I don't see that.  I see response amplitude vs omega.  Sorry if that sounds dense but I don't know what you're doing.

But in this case the black points are given are 'a,b,c,d' have to be extracted.

One of the major challenges is that is that there are multiples "y" for one "x" so the function cannot associate uniquely the two set of data.

Moreover it is not possible to have an explicit 'y' with respect to 'x'.

OK.  Again, that makes it seem to me like it is not a data fitting problem at all.  You have x and y and want to find a, b, c, and d such that your expression is satisfied.  That's more like a root-finding problem.

In some cases the lmfit function works and it can extract parameters correctly but, in case for two consequently x (x2-x1) there is a big gap (y2-y1) The model explode and gives a 'd' value too big.

An example is in the picture 'lmfit_badresult'

I don't fully understand the plot (probably similar to above -- are x and y related to amplitude and frequency?), but it looks like it must be a very bad set of initial values.

Why are you re-doing the fit?  What is that trying to accomplish?

I use a second fitting to confront the first fitting with the second to make sure that the algorithm reached the steady state and there is no a significant variation respect to the initial guess.

I think that won't actually do much.  The solution found is actually deterministic and generally includes no randomness, at least for most methods and you didn't specify which fitting method you were using.

Hope that helps,

--Matt