Completely wrong result when I use weighting function

[Shankar Dutt]

Hi everyone,

Thanks a lot in advance for your help. I have tried a lot before posting here but I am unable to understand. I have experimental data with q, r and s (q-x axis values, r - y-axis values and s is the errors). I have attached the data as well. If I try to fit without the weighting factor(s), I am getting what I expect but with weighting factor added as shown, I am getting completely wrong results. Don't know what am I doing wrong?

import matplotlib.pyplot as plt

import numpy as np

from lmfit import Minimizer, Parameters, fit_report, minimize,  Model, report_fit

import models.Models1D as mdls

import scipy.special as spe

from numpy import exp, sqrt, pi

from scipy.integrate import quad_vec

filepath = 'IGOR_1168.txt'

data = np.loadtxt(filepath)

q = data[:, 0]

r = data[:, 1]

s = data[:, 2]

p = Parameters()

p.add('Amplitude', value=0.1, vary=True)

p.add('Radius', value=600, vary=True)

p.add('Polydispersity', value=0.03, vary=False)

p.add('Background', value=0, vary=False)

p.add('Roughness', value=0, vary=False)

def hc(amp,radius,pd,background,roughness):

    s = pd*radius

    va = radius-4*s

    vb = radius+4*s   

    

    def dist_fun(z):

        return (1 / (sqrt(2 * pi) * s)) * exp(-0.5*(((z - radius) / s)**2))

    dist_int = quad_vec(dist_fun, va, vb)[0]

    def fun(z):

        return ((1/q) * (z * spe.j1(q * z)))**2 * dist_fun(z)

    def integral_hc():

        return quad_vec(fun,va,vb)[0]

    

    return (integral_hc() / dist_int) * amp ** 2 * np.exp(-(roughness*q)**2)  + background

def fcn2mmin(p, x, data = None, eps = None):

    p = p.valuesdict()

    amp = p['Amplitude']

    radius = p['Radius']

    pd = p['Polydispersity']

    background = p['Background']

    roughness = p['Roughness']

    model = hc(amp,radius,pd,background,roughness)

    if data is None:

        return model

    return (model-data)/eps

x = q

data = r

eps = s

# do fit

minner = Minimizer(fcn2mmin, params=p, fcn_args=(x,), fcn_kws={'data':data,'eps':eps})

result = minner.leastsq()

# calculate final result

final = data + result.residual

print(result.params)

# write error report

fr = fit_report(result)

print(fr)

font = {'family': 'serif',

            'color': 'darkred',

            'weight': 'normal',

            'size': 16,

            }

plt.loglog(x, data,  label='DATA',linewidth=0.7)

plt.loglog(x, final, 'r', label='FIT',linewidth=1.5)

plt.xlabel(r'q', fontdict=font)

plt.ylabel(r'I', fontdict=font)

plt.title('Hard Cylinder Model', fontdict=font)

plt.legend()

plt.show()

The code using the model class gives expected results. Is it my mistake? or some bug in lmfit? The code is as follows:

 

def mdl(x, Amplitude, Background, Radius, Polydispersity, Roughness):
     return mdls.hc_r_model(x, Amplitude, Background, Radius, Polydispersity, Roughness, wavelength)
 
 fmodel = Model(mdl)
 params = fmodel.make_params()
 params = p
 result = fmodel.fit(r, params, x=q, weights = s)
 fr = result.fit_report()
 print(fr)
 final = result.best_fit

[Matt Newville]

On Fri, Apr 24, 2020 at 6:42 PM Shankar Dutt <shankar...@gmail.com> wrote:

Hi everyone,

Thanks a lot in advance for your help. I have tried a lot before posting here but I am unable to understand. I have experimental data with q, r and s (q-x axis values, r - y-axis values and s is the errors). I have attached the data as well. If I try to fit without the weighting factor(s), I am getting what I expect but with weighting factor added as shown, I am getting completely wrong results. Don't know what am I doing wrong?

It's hard to know without knowing what you mean by "completely wrong results".   Like, did the fit run to completion and give results you think are not possible, or did the fit not run at all?  You are dividing by 'eps', so an obvious question would be: are any of these values 0?  That would definitely be a problem.

I haven't looked at your data, but I see that you would be plotting on a log-log scale.  Just as a general rule-of-thumb: if you are plotting on a log-log or semi-log scale, you might find it useful to fit log(data) to log(model).  I don't know if that is what is causing the trouble you are seeing,

--Matt

The code without the weighting function is 

def fcn2mmin(p, x, data = None):

    p = p.valuesdict()

    amp = p['Amplitude']

    radius = p['Radius']

    pd = p['Polydispersity']

    background = p['Background']

    roughness = p['Roughness']

    model = hc(amp,radius,pd,background,roughness)

    if data is None:

        return model

    return (model-data)

x = q

data = r

eps = s

# do fit

minner = Minimizer(fcn2mmin, params=p, fcn_args=(x,), fcn_kws={'data':data})

result = minner.leastsq()

# calculate final result

final = data + result.residual

print(result.params)

# write error report

fr = fit_report(result)

print(fr)

font = {'family': 'serif',

            'color': 'darkred',

            'weight': 'normal',

            'size': 16,

            }

plt.loglog(x, data,  label='DATA',linewidth=0.7)

plt.loglog(x, final, 'r', label='FIT',linewidth=1.5)

plt.xlabel(r'q', fontdict=font)

plt.ylabel(r'I', fontdict=font)

plt.title('Hard Cylinder Model', fontdict=font)

plt.legend()

plt.show()

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/503f6700-f1bb-4f9e-b15e-c2051c033246%40googlegroups.com.

[Shankar Dutt]

Dear Matt,

Thanks a lot for your reply. Sorry for late reply. Had lost access due to the internet and due to coronavirus situation, took sometime to fix it. Anyhow, by completely wrong result, I mean it just random parameter results and minimize.residue also final = data + result.residual gave me wrong result. Using model class, I am getting what I expected.  Also, eps has no zero values.  I can attach graphs, if u think that will help. Also, can you please clarify the point you mentioned earlier " if you are plotting on a log-log or semi-log scale, you might find it useful to fit log(data) to log(model)". Why so? and what is the most efficient way to do so?


To unsubscribe from this group and stop receiving emails from it, send an email to lmfi...@googlegroups.com.

[Matt Newville]

On Tue, Apr 28, 2020 at 2:59 PM Shankar Dutt <shankar...@gmail.com> wrote:

Dear Matt,

Thanks a lot for your reply. Sorry for late reply. Had lost access due to the internet and due to coronavirus situation, took sometime to fix it. Anyhow, by completely wrong result, I mean it just random parameter results and minimize.residue also final = data + result.residual gave me wrong result. Using model class, I am getting what I expected. I can attach graphs, if u think that will help.

You did not answer my earlier questions.

 

Also, can you please clarify the point you mentioned earlier " if you are plotting on a log-log or semi-log scale, you might find it useful to fit log(data) to log(model)". Why so? and what is the most efficient way to do so?

Well, if you're plotting log(ydata) vs. xdata or vs. log(xdata), then you might want your model function to return log(model) and fit log(ydata) to that.  

That will give greater emphasis to the smaller values.  

--Matt