#!/usr/bin/env python
# coding: utf-8

# # iPython Cookbook - Curve Fitting

# In[1]:


import numpy as np
import scipy as sp
import matplotlib.pyplot as plt


# ## Generating some data to play with
# We first generate some data to play with. So in the first step we choose a functional relationship $func_0$ between our $x$ and $y$ that is parameterised by three parameters $a,b,c$...

# In[2]:


a,b,c=(1,2,1)

def func0 (x,a,b,c):
    return a*exp(-b*x)+c


# ...we now choose the area over which we want to plot the function $x_{min}\ldots x_{max}$ and the number of points we want to generate $N$ and generate the undisturbed curve $yvals_0$ over a random selection of $x$ coordinates $xvals$.

# In[3]:


xmin,xmax = (0,1)
N = 500
xvals = np.random.uniform(xmin, xmax, N)
yvals0 = func0(xvals,a,b,c)
plt.plot(xvals, yvals0, '+')
plt.show()


# In a second step we generate the error term $err$ and add it to the previously computed values, storing the results in $yvals$

# In[4]:


sig = 0.05
err = sig * np.random.standard_normal(N)
yvals = yvals0 + err
plt.plot(xvals, yvals, '+')
plt.show()


# ##Fitting the curve
# We now import the relvant module, and we now choose the functional relationship $func$ between our $x$ and $y$ that is parameterised by three parameters $a,b,c$ of the curve we want to fit. We here use the same relationship as above, but obviously we can fit any curve we want (eg, linear, quadratic etc). We also have to provide an initial guess for the paramters, here $a_i,b_i,c_i$.

# In[5]:


# import the curve fitting module and standard imports
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit 

# choose the function to be fitted... 
def func (x,a,b,c):
    return a*exp(-b*x)+c

# ...and provide initial estimates for the parameters
a0,b0,c0 = (0.5,0.5,0.5)


# We now execute the curve fit and plot the results. The coefficients will be in the tuple $coeffs=(a,b,c)$ and $fiterr$ will contain an error estimate. If $fiterr$ is NaN (or the algorithm terminates with too many iterations) this means that the algorithm did not converge, which means we need to adapt the initial values (see below)

# In[6]:


# exectute the curve fit...
coeffs, fiterr = curve_fit(func, xvals, yvals, p0=(a0,b0,c0))

# ...and plot the results
print ("a=%s, b=%s, c=%s" % (coeffs[0], coeffs[1], coeffs[2]))
plt.plot(xvals,yvals, '+')
plt.plot(xvals,func(xvals,*coeffs),'r.')
plt.show()


# In order to find the initial values we might want to plot the function with some initial parameters $a_t,b_t,c_t$ to allow for a rough manual fit

# In[7]:


# manually fit the curve to obtain a viable set of starting parameters
at,bt,ct = (0.5,0.5,0.5)
plt.plot(xvals,yvals, '+')
plot(xvals,func(xvals,at,bt,ct), 'g.')
print ("a=%s, b=%s, c=%s" % (at,bt,ct))


# ## Licence and version
# *(c) Stefan Loesch / oditorium 2014; all rights reserved 
# ([license](https://github.com/oditorium/blog/blob/master/LICENSE))*

# In[8]:


import sys
print(sys.version)