chi-squared calculation different than what I expected

[jprince58]

I'm using lmfit to do a curve-fitting to a PDE solved using a finite-difference method.

When i take the result from my fitting using result=lmfit.minimize(), the chi-squared value (pulled from result.chisqr), is different than what I am calculating if I do the chi-squared calculation myself.

Not sure what is going on. Can link the github repo containing the scripts as necessary.

[newville]

@jprince58 Um:

Not sure what is going on.

Imagine how we feel. Some random person on the internet literally said "I did a calculation and it is different from what I expected". Good grief.

Can link the github repo containing the scripts as necessary.

Well, yeah. If you don't tell us what you did, we do not what you did.

[jprince58]

Fair. I'm not an expert on these kinds of forums, so I wasn't sure how to start the thread.

This repo has all the associated files.

https://github.com/jprince58/Bioaccumulation_Nanoparticles/tree/Quantitative_fit/N12

I wish I could just tell you, "run this code and you'll replicate the issue", but the scripts are set-up to run on only my machine, so you would need to change the directories referenced. The first thing that needs to be changed are lines 19 and 20 of "N12-Fit_Code", and should reference any directory one wants on the machine running the script, as long as it can hold files. The "direct_export_path" variable is where a document summarizing the most recent code execution lives. Additionally, lines 56 and 126 of the same script would need to be changed to reference the attached .csv file. Additonally, whatever directory the script is going to be run in will need a plain .txt file added named "counter_file_1", with a single number entry. This corresponds with lines 48-53 of "N12-Fit_Code".

I'll just focus on the part of the code with the issue. It's the "N12-Fit_Code" script in the directory referenced. Lines 94-100 look as such:

Starting at line 94

fit_results = minimize(residual_calc, params, method='lbfgsb')  # fitting for model   
RSS=fit_results.chisqr 
p=9 #Number of parameters
N=227 #Number of data-points used in fit (looked at excel sheet)
AICc=2*p+N*(np.log(2*3.14*RSS/N)+1)

I suspect that the RSS which is calculated in line 95 of "N12-Fit_Code" is off because in the file "N12_exp_data_fitter", the RSS was calculated independently in line 94 (calculated with the same values as those passed along in "residual_calc" used in the minimize function)

This can be seen most easily if N12-Fit_Code is ran with the following parameter set:

params.add('h', value=0.05/2, vary=False)
params.add('tol', value=10**(-8), vary=False)
params.add('t1', value=0, vary=False)
params.add('t2', value=15, vary=False)
params.add('nx', value=50, vary=False)
# params.add('omega', value=0.8, min=0.5, max=1.)
params.add('omega', value=1., vary=False)
# params.add('mu', value=0.01, min=0., max=10.)
params.add('mu', value=0.01, vary=False)
# params.add('nu', value=10, min=0.01, max=100.)
params.add('nu', value=25, vary=False)
# params.add('phi_min', value=0.3, min=0.1, max=1)
params.add('phi_min', value=0.25, vary=False)
# params.add('phi_max', value=0.8, min=0.1, max=1)
params.add('phi_max', value=0.8, vary=False)
# params.add('kappa', value=20, min=0.1, max=50.)
params.add('kappa', value=20, vary=False)
params.add('a', value=2.5, min=1., max=50.) 
# params.add('a', value=1.05, vary=False)  
# params.add('b', value=25, min=0, max=100.)
params.add('b', value=90, vary=False)
# params.add('c', value=5.3, min=0.01, max=10.)
params.add('c', value=5.3, vary=False)
# params.add('Kp', value=5, min=0.01, max=10.)
params.add('Kp', value=1, vary=False)
params.add('kconv', value=188.8, vary=False)

If you run the code with this parameter set, and you find the exported .docx file in the directory named in line 20, you should hopefully see something like this. The RSS typed out in the doc is that calculated using lmfit. The one in the figure itself is the one calculated in "N12_exp_data_fitter" line 94. Obviously they are different.

Attached is the curve (more like set of curves) trying to be fit.

tseng_fits_Fig2B_Cy5cipro_incubation_bump.csv

I understand this probably isn't enough information. Please let me know any other information needed.

[reneeotten]

very unlikely that anybody will actually go through all of this.

Show us the code where you calculate the chisqr yourself and how different the value is from what lmfit gives you.

[jprince58]

Calculated myself is in N12_exp_data_fitter, lines 93 and 94:

        res=collated_results[:,2]-collated_results[:,3]
        RSS=np.nansum(res**2)

The difference in values can be seen in the attached figure above.

[newville]

If you cannot make a minimal working example, we cannot help you. See https://stackoverflow.com/help/minimal-reproducible-example.