Skip to content
Snippets Groups Projects
Commit 19b57b94 authored by Saad Jbabdi's avatar Saad Jbabdi
Browse files

cleaned

parent b99461f1
No related branches found
No related tags found
1 merge request!21Saad edits
Expand all Hide whitespace changes
Inline Side-by-side
talks/matlab_vs_python/rbf/.ipynb_checkpoints/fit_model-checkpoint.ipynb deleted 100644 → 0
+ 0
− 157
View file @ b99461f1
%% Cell type:code id: tags:
``` python
# Fit a model to some data
# Model is:
# prediction = M0 * exp(-TE/T2)*(1-exp(-TR/T1))
# where M0,T1,T2 are unknown parameters and TE/TR are experimental parameters
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
TEs = np.array([10,40,60,80]) # TE values in ms
TRs = np.array([.5,1,1.5,2]) # TR in seconds
# All combinations of TEs/TRs
combinations = np.array([(x,y) for x in TEs for y in TRs])
TEs,TRs = combinations[:,0],combinations[:,1]
# function for our model
def forward(p):
M0,T1,T2 = p
return M0*np.exp(-TEs/T2)*(1-np.exp(-TRs/T1))
# simulate data using model
true_p = [100,.8,50]
data = forward(true_p)
data = data + np.random.randn(data.size)
```
%% Output
[<matplotlib.lines.Line2D at 0x121ca4278>]
%% Cell type:code id: tags:
``` python
# Now for the fitting
# we need a cost function:
def cf(p):
pred = forward(p)
return np.mean((pred-data)**2)/2
# always a good idea to calculate gradient
def forward_deriv(p):
M0,T1,T2 = p
E1,E2 = np.exp(-TEs/T2),np.exp(-TRs/T1)
dfdM0 = E2*(1-E1)
dfdT1 = M0*E2*(-E1/T1**2)
dfdT2 = M0*(E2/T2**2)*(1-E1)
return np.array([dfdM0,dfdT1,dfdT2])
def gradient(p):
pred = forward(p)
deriv = forward_deriv(p)
return np.mean( deriv * (pred-data)[None,:],axis=1)
# get ready to minimize
p0 = [200,70,1000] # some random guess
method = 'TNC'
arguments = {'x0':p0,'method':method,'jac':gradient}
result = minimize(cf,**arguments)
```
%% Cell type:code id: tags:
``` python
plt.plot(data)
plt.plot(forward(result.x))
```
%% Output
[<matplotlib.lines.Line2D at 0x121f34eb8>]
%% Cell type:code id: tags:
``` python
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment