Using a guess with scipy curve_fit

Question

I have a function that I want to curve fit with knowing the error of the curve fit. I'm trying to use scipy.optimize.curve_fit to do this but am running into problem. Right now my code is:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

pi = np.pi
sqrt = np.sqrt
mean = np.mean

A = 1
T_2 = 100
nu_0 = 10
phi_0 = 0
n = .001
nu_s = 500
T = 1000

t = np.linspace(0, T, num = (nu_s*T))

def S_n(A,t,T_2,nu_0,phi_0,n,nu_s,T):
    return (A/np.sqrt(2))*np.exp(-t/T_2)*np.cos(2*pi*nu_0*t+phi_0)

S = S_n(A,t,T_2,nu_0,phi_0,n,nu_s,T) + np.random.normal(0, n, nu_s*T)

guess = np.array([A,T_2,nu_0,phi_0,n,nu_s,T])
print guess

popt, pcov = curve_fit(S_n,t,S, guess)
print popt
perr = sqrt(np.diag(pcov))
print perr

Which gives me infinite error. I'm not sure I am doing my guess correctly because in the equation, everything remains constant except t, so do I leave t out of my guess because it is no longer an array because t is a sequence. When I leave t out of the guess, which I am doing here I receive values for each variable to be far from the values I give initially with infinite error. If I include t in the guess, then I get an error.

Answer 1

You didn't take the order of the parameters to curve_fit into account:

Definition: curve_fit(f, xdata, ydata, p0=None, sigma=None, **kw)

Docstring: Use non-linear least squares to fit a function, f, to data.

Assumes ydata = f(xdata, *params) + eps

Parameters

f : callable The model function, f(x, ...). It must take the independent variable as the first argument and the parameters to fit as separate remaining arguments.

If you make a function:

def Sm(t, A,T_2,nu_0,phi_0,n,nu_s,T):
    return S_n(A, t, T_2,nu_0,phi_0,n,nu_s,T)

(notice the order of the first 2 parameters has changed) and pass that to curve_fit, it will work.

It is probably more pythonic to clean your original function though:

def S_n(t, A, T_2, nu_0, phi_0, n, nu_s, T):
    return (A/np.sqrt(2))*np.exp(-t/T_2)*np.cos(2*pi*nu_0*t+phi_0)

S = S_n(t, A, T_2, nu_0, phi_0, n, nu_s, T) + np.random.normal(0, n, nu_s*T)

Then you can pass S_n to curve_fit unchanged, as well as S and guess.

To be clear, the following code produces the desired fit:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def S_n(t, amplitude, sigma,freq0,phase):
    return (amplitude/np.sqrt(2))*np.exp(-t/sigma)*np.cos(2*np.pi*freq0*t+ phase)

amplitude = 1
sigma = 100
freq0 = 10
phase = 0
n = .001
sample_freq = 500
period = 1000

t = np.linspace(0, period, num = (sample_freq*period))
S = S_n(t, amplitude, sigma,freq0,phi_0) + np.random.normal(0, n, sample_freq*period)

guess = np.array([amplitude, sigma, freq0, phase ])
print guess
popt, pcov = curve_fit(S_n,t,S, guess)
print popt
print(np.all(np.isfinite(pcov)))
# output
[  1 100  10   0]
[  1.00000532e+00   1.00000409e+02   1.00000000e+01  -1.49076430e-05]
True