In [13]:

    

# Load Modules
import numpy as np
from scipy.optimize import curve_fit
%matplotlib inline  
import matplotlib.pyplot as plt


    

In [14]:

    

# Define function to be used
def T2decay(xdata, T1,Mz):
    R1=1.0/T1
    return Mz * np.exp(-xdata*R1)

def plot(xdata,ydata1,ydata2):
    plt.plot(xdata,ydata1)
    plt.scatter(xdata,ydata2,color='r')
    plt.show()


    

In [15]:

    

# Simulate data
NPoints=20
TR=np.linspace(0,0.2,NPoints)
T2=.050
Mz=1.0
Signal=T2decay(TR, T2,Mz)


    

In [16]:

    

# Add Gaussian noise
mu=0.0; sigma=0.1;
Signal_with_noise= Signal + np.random.normal(mu, sigma, NPoints)


    

In [17]:

    

# Plot
plot(TR,Signal,Signal_with_noise)


    

In [18]:

    

# Perform Non-Linear curve fitting with bounds
# using SciPy
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html
initial_guess=[1,.1]
pars_hat, cov = curve_fit(T2decay, TR, Signal_with_noise,p0=initial_guess,bounds=(0, 1.0))


    

In [19]:

    

# estimate Yhat
Yhat = T2decay(TR, pars_hat[0],pars_hat[1])


    

In [20]:

    

# Plot estimated data
plot(TR,Signal,Yhat)


    

In [21]:

    

# Remember to Install if you don't have it
from lmfit import minimize, Parameters


    

In [22]:

    

# LMFIT requires a function that returns the objective function. In this case the sum of the square of the residuals
def residual(params, xdata, ydata_observed):
    R2 = 1.0/params['T2']
    Mz = params['Mz']
    model = Mz * np.exp(-xdata*R2)
    return (ydata_observed-model)

params = Parameters()
params.add('T2', value=0.1,min=0.0,max=1.0)
params.add('Mz', value=1.0,min=0.0,max=1.0)


    

In [40]:

    

# run fitting
lmfit = minimize(residual, params, args=(TR, Signal_with_noise))


    

In [49]:

    

# estimate Yhat
T2pred=lmfit.params['T2'].value
Mzpred=lmfit.params['Mz'].value
Yhat2 = T2decay(TR, T2pred,Mzpred)


    

In [50]:

    

# Plot estimated data
plot(TR,Signal,Yhat2)


    

In [ ]: