Why does scipy.optimize.curve_fit not fit correctly to the data?

Question

I've been trying to fit a function to some data for a while using scipy.optimize.curve_fit but I have real difficulty. I really can't see any reason why this wouldn't work.

# encoding: utf-8
from __future__ import (print_function,
                        division,
                        unicode_literals,
                        absolute_import,
                        with_statement)
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as mpl

x, y, e_y = np.loadtxt('data.txt', unpack=True)

def f(x, a, k):
    return (1/(np.sqrt(1 + a*((k-x)**2))))

popt, pcov = curve_fit(f, x, y, maxfev = 100000000)

mpl.plot(x, f(x, *popt), 'r-', label='Fit')
mpl.plot(x, y, 'rx', label='Original')
mpl.legend(loc='best')
mpl.savefig('curve.pdf')
print(popt)

# correct values which should be calculated
# a=0.003097
# k=35.4

Here is the plot-image which is produced by upper code:

data.txt:
#x      y       e_y
4.4     0.79    0.13
19.7    4.9     0.8
23.5    7.3     1.2
29.7    17      2.79
30.7    21.5    3.52
34      81      13.28
34.6    145     23.77
35.4    610     100
36.3    115     18.85
38.1    38      6.23
43.7    14      2.3
56.2    6.2     1.02
64.7    4.6     0.75
79.9    3.2     0.52
210     0.98    0.16

Answer 1

Firstly try not to increase maxfev so large, this is usually a sign something else is going wrong! Playing around I can get a fit by the following addition:

def f(x, b, a, k):
    return (b/(np.sqrt(1 + a*((k-x)**2))))

popt, pcov = curve_fit(f, x, y, p0=[20, 600.0, 35.0])

Firstly give the fitting function you have given has a maximum of 1, since the peak in your data is 600, it will never fit. So I added an overall factor b. Secondly , try to help poor old curve_fit out. If by eye you can see it peaks at x~35 then tell it through the p0. This requires some intuition as to how the function works but is very important if your going to use a curve fitting function.