3

Running the following code,

x = np.array([50.849937, 53.849937, 56.849937, 59.849937, 62.849937, 65.849937, 68.849937, 71.849937, 74.849937, 77.849937, 80.849937, 83.849937, 86.849937, 89.849937, 92.849937])
y = np.array([410.67800, 402.63800, 402.63800, 386.55800, 330.27600, 217.71400, 72.98990, 16.70860, 8.66833, 40.82920, 241.83400, 386.55800, 394.59800, 394.59800, 402.63800])
def f(om, a, i , c):
       return a - i*np.exp(- c* (om-74.)**2)
par, cov = curve_fit(f, x, y)
stdev = np.sqrt(np.diag(cov) )

produces this Graph,

enter image description here

With the following parameters and standard deviation:

par =   [ 4.09652163e+02, 4.33961227e+02, 1.58719772e-02]
stdev = [ 1.46309578e+01, 2.44878171e+01, 2.40474753e-03]

However, when trying to fit this data to the following function:

def f(om, a, i , c, omo):
       return a - i*np.exp(- c* (om-omo)**2)

It doesn't work, it produces a standard deviation of

stdev = [inf, inf, inf, inf, inf]

Is there any way to fix this?

Ed Smith
  • 12,716
  • 2
  • 43
  • 55
iqopi
  • 141
  • 1
  • 2

2 Answers2

3

It looks like it isn't converging (see this and this). Try adding an initial condition,

par, cov = curve_fit(f, x, y, p0=[1.,1.,1.,74.])

which results in the

par = [ 4.11892318e+02, 4.36953868e+02, 1.55741131e-02, 7.32560690e+01])
stdev = [ 1.17579445e+01, 1.94401006e+01, 1.86709423e-03, 2.62952690e-01]
Community
  • 1
  • 1
Ed Smith
  • 12,716
  • 2
  • 43
  • 55
  • had to set all condition: 'startfirst = [86.5, -438., -411., 73.5]', didn't work with just one, but thanks – iqopi Sep 27 '16 at 08:11
2

You can calculate the initial condition from data:

%matplotlib inline
import pylab as pl

import numpy as np
from scipy.optimize import curve_fit

x = np.array([50.849937, 53.849937, 56.849937, 59.849937, 62.849937, 65.849937, 68.849937, 71.849937, 74.849937, 77.849937, 80.849937, 83.849937, 86.849937, 89.849937, 92.849937])
y = np.array([410.67800, 402.63800, 402.63800, 386.55800, 330.27600, 217.71400, 72.98990, 16.70860, 8.66833, 40.82920, 241.83400, 386.55800, 394.59800, 394.59800, 402.63800])

def f(om, a, i , c, omo):
       return a - i*np.exp(- c* (om-omo)**2)

par, cov = curve_fit(f, x, y, p0=[y.max(), y.ptp(), 1, x[np.argmin(y)]])
stdev = np.sqrt(np.diag(cov) )

pl.plot(x, y, "o")
x2 = np.linspace(x.min(), x.max(), 100)
pl.plot(x2, f(x2, *par))
HYRY
  • 94,853
  • 25
  • 187
  • 187