In which we explore temporal derivatives for dealing with timing issues.
In [1]:
%matplotlib inline
In [2]:
import numpy as np
import matplotlib.pyplot as plt
Get HRF models:
In [3]:
from nipy.modalities.fmri import hrf
In [4]:
# Times from 0 to 30 seconds at 0.1 second intervals
time_30s = np.arange(0, 30, 0.1)
# HRF for these times
canonical_hrf = hrf.spmt(time_30s)
# SPM's own temporal derivative
canonical_dhrf = hrf.dspmt(time_30s)
# Something less unexpected for the temporal derivative
# canonical_dhrf = np.gradient(canonical_hrf)
In [5]:
plt.plot(time_30s, canonical_hrf)
plt.plot(time_30s, canonical_dhrf)
Out[5]:
Make lots of simulated HRF shapes, offset by 1 time delta
In [6]:
time_60s = np.arange(0, 60, 0.1)
long_hrf = hrf.spmt(time_60s)
In [7]:
import scipy.linalg as spl
Toeplitz routine is a trick for making HRFs offset at 0 seconds, 0.1 seconds, 0.2 seconds ...
In [8]:
hrf_collection = spl.toeplitz(long_hrf, np.zeros_like(time_30s))
In [9]:
plt.imshow(hrf_collection)
Out[9]:
In [10]:
time_15s = time_30s[time_30s < 15]
n_time_15s = len(time_15s)
n_time_15s
Out[10]:
In [11]:
time_pm15s = time_30s - 15
time_pm15s[n_time_15s]
Out[11]:
In [12]:
hrf_model = hrf_collection[:, n_time_15s]
dhrf_model = hrf_model.copy()
dhrf_model[n_time_15s:n_time_15s+len(canonical_dhrf)] = canonical_dhrf
In [13]:
plt.plot(time_60s, hrf_model)
plt.plot(time_60s, dhrf_model)
Out[13]:
In [14]:
np.all(hrf_model == hrf_collection[:, n_time_15s])
Out[14]:
Make the design matrix
In [15]:
X = np.vstack((hrf_model, dhrf_model, np.ones_like(hrf_model))).T
In [16]:
plt.imshow(X, aspect='auto', interpolation='nearest', cmap='gray')
Out[16]:
Fit the offset data with this design:
In [17]:
Y = hrf_collection
B = spl.pinv(X).dot(hrf_collection)
fitted = X.dot(B)
E = Y - fitted
In [18]:
rmse = np.sqrt((E ** 2).sum(axis=0))
rmse.shape
Out[18]:
In [19]:
plt.plot(time_pm15s, rmse)
Out[19]:
In [20]:
time_pm4s_i = (time_pm15s >= -4) & (time_pm15s < 4)
plt.plot(time_pm15s[time_pm4s_i], rmse[time_pm4s_i])
Out[20]:
In [21]:
tm1_i = np.nonzero(time_pm15s == -1)[0]
t0_i = np.nonzero(time_pm15s == 0)[0]
tp5_i = np.nonzero(time_pm15s == 0.5)[0]
t1_i = np.nonzero(time_pm15s == 1)[0]
t2_i = np.nonzero(time_pm15s == 2)[0]
tm1_i, t0_i, tp5_i, t1_i, t2_i
Out[21]:
In [22]:
def plot_at_index(i):
beta_hrf, beta_dhrf, beta_mean = B[:, i]
plt.plot(time_60s, hrf_collection[:, i], label='y')
plt.plot(time_60s, hrf_model * beta_hrf, label='fitted hrf')
plt.plot(time_60s, dhrf_model * beta_dhrf, label='fitted dhrf')
plt.plot(time_60s, fitted[:, i], label='fitted sum')
plt.plot(time_60s, E[:, i], label='error')
plt.legend()
plot_at_index(t0_i)
In [23]:
plot_at_index(tp5_i)
In [24]:
# For SPM temporal derivative hrf_t0 + -(hrft0 - hrft1)
plot_at_index(t1_i)
In [25]:
plot_at_index(tm1_i)
In [26]:
plot_at_index(t2_i)
In [27]:
# Parameter for the first (HRF) regressor
plt.plot(time_pm15s, B[0,:])
# Parameter for the second (temporal derivative) regressor
plt.plot(time_pm15s, B[1,:])
Out[27]:
How about orthogonality?
In [28]:
canonical_hrf.dot(canonical_dhrf)
Out[28]:
In [29]:
def cos_angle(v1, v2):
""" Cosine of angle between vectors
0 means orthogonal, 1 means co-linear
"""
len_v1 = np.sqrt(v1.dot(v1))
len_v2 = np.sqrt(v2.dot(v2))
return v1.dot(v2) / (len_v1 * len_v2)
cos_angle(canonical_hrf, canonical_dhrf)
Out[29]:
Not quite. But now let's do some convolution.
In [30]:
neural_signal = np.zeros_like(time_60s)
neural_signal[time_60s < 10] = 1
plt.plot(time_60s, neural_signal)
Out[30]:
In [31]:
n_60s = len(time_60s)
hemo_regressor_hrf = np.convolve(canonical_hrf, neural_signal)[:n_60s]
hemo_regressor_dhrf = np.convolve(canonical_dhrf, neural_signal)[:n_60s]
len(hemo_regressor_hrf)
Out[31]:
In [32]:
plt.plot(time_60s, hemo_regressor_hrf)
plt.plot(time_60s, hemo_regressor_dhrf)
Out[32]:
Are the regressors orthogonal?
In [33]:
np.dot(hemo_regressor_hrf, hemo_regressor_dhrf)
Out[33]:
In [34]:
cos_angle(hemo_regressor_hrf, hemo_regressor_dhrf)
Out[34]:
How about longer block?
In [35]:
neural_signal[time_60s < 30] = 1
hemo_regressor_hrf = np.convolve(canonical_hrf, neural_signal)[:n_60s]
hemo_regressor_dhrf = np.convolve(canonical_dhrf, neural_signal)[:n_60s]
plt.plot(time_60s, hemo_regressor_hrf)
plt.plot(time_60s, hemo_regressor_dhrf)
Out[35]:
In [36]:
cos_angle(hemo_regressor_hrf, hemo_regressor_dhrf)
Out[36]:
Shorter?
In [37]:
neural_signal[:] = 0
neural_signal[time_60s < 1] = 1
hemo_regressor_hrf = np.convolve(canonical_hrf, neural_signal)[:n_60s]
hemo_regressor_dhrf = np.convolve(canonical_dhrf, neural_signal)[:n_60s]
plt.plot(time_60s, hemo_regressor_hrf)
plt.plot(time_60s, hemo_regressor_dhrf)
Out[37]:
In [38]:
cos_angle(hemo_regressor_hrf, hemo_regressor_dhrf)
Out[38]:
In [39]:
def cos_angles_for_t(t, times=time_60s):
neural_signal = np.zeros_like(times)
n = len(neural_signal)
neural_signal[times < t] = 1
hemo_regressor_hrf = np.convolve(canonical_hrf, neural_signal)[:n]
hemo_regressor_dhrf = np.convolve(canonical_dhrf, neural_signal)[:n]
return cos_angle(hemo_regressor_hrf, hemo_regressor_dhrf)
In [40]:
cos_angles = []
for t in time_30s:
cos_angles.append(cos_angles_for_t(t))
In [41]:
plt.plot(time_30s, cos_angles)
Out[41]:
