In which we explore temporal derivatives for dealing with timing issues.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
Get HRF models:
from nipy.modalities.fmri import hrf
# 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)
plt.plot(time_30s, canonical_hrf)
plt.plot(time_30s, canonical_dhrf)
[<matplotlib.lines.Line2D at 0x108ad8c90>]
Make lots of simulated HRF shapes, offset by 1 time delta
time_60s = np.arange(0, 60, 0.1)
long_hrf = hrf.spmt(time_60s)
import scipy.linalg as spl
Toeplitz routine is a trick for making HRFs offset at 0 seconds, 0.1 seconds, 0.2 seconds ...
hrf_collection = spl.toeplitz(long_hrf, np.zeros_like(time_30s))
plt.imshow(hrf_collection)
<matplotlib.image.AxesImage at 0x108d24f50>
time_15s = time_30s[time_30s < 15]
n_time_15s = len(time_15s)
n_time_15s
150
time_pm15s = time_30s - 15
time_pm15s[n_time_15s]
0.0
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
plt.plot(time_60s, hrf_model)
plt.plot(time_60s, dhrf_model)
[<matplotlib.lines.Line2D at 0x108cec550>]
np.all(hrf_model == hrf_collection[:, n_time_15s])
True
Make the design matrix
X = np.vstack((hrf_model, dhrf_model, np.ones_like(hrf_model))).T
plt.imshow(X, aspect='auto', interpolation='nearest', cmap='gray')
<matplotlib.image.AxesImage at 0x10a1f3710>
Fit the offset data with this design:
Y = hrf_collection
B = spl.pinv(X).dot(hrf_collection)
fitted = X.dot(B)
E = Y - fitted
rmse = np.sqrt((E ** 2).sum(axis=0))
rmse.shape
(300,)
plt.plot(time_pm15s, rmse)
[<matplotlib.lines.Line2D at 0x10b390590>]
time_pm4s_i = (time_pm15s >= -4) & (time_pm15s < 4)
plt.plot(time_pm15s[time_pm4s_i], rmse[time_pm4s_i])
[<matplotlib.lines.Line2D at 0x10b49b490>]
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
(array([140]), array([150]), array([155]), array([160]), array([170]))
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)
plot_at_index(tp5_i)
# For SPM temporal derivative hrf_t0 + -(hrft0 - hrft1)
plot_at_index(t1_i)
plot_at_index(tm1_i)
plot_at_index(t2_i)
# 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,:])
[<matplotlib.lines.Line2D at 0x10c1acf90>]
How about orthogonality?
canonical_hrf.dot(canonical_dhrf)
0.099261698796546796
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)
0.16767552736516803
Not quite. But now let's do some convolution.
neural_signal = np.zeros_like(time_60s)
neural_signal[time_60s < 10] = 1
plt.plot(time_60s, neural_signal)
[<matplotlib.lines.Line2D at 0x10c640750>]
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)
600
plt.plot(time_60s, hemo_regressor_hrf)
plt.plot(time_60s, hemo_regressor_dhrf)
[<matplotlib.lines.Line2D at 0x10c678890>]
Are the regressors orthogonal?
np.dot(hemo_regressor_hrf, hemo_regressor_dhrf)
191.58719003912032
cos_angle(hemo_regressor_hrf, hemo_regressor_dhrf)
0.09848538644886784
How about longer block?
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)
[<matplotlib.lines.Line2D at 0x10c694890>]
cos_angle(hemo_regressor_hrf, hemo_regressor_dhrf)
0.052984565515914336
Shorter?
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)
[<matplotlib.lines.Line2D at 0x10d05ab10>]
cos_angle(hemo_regressor_hrf, hemo_regressor_dhrf)
0.16580167893517508
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)
cos_angles = []
for t in time_30s:
cos_angles.append(cos_angles_for_t(t))
-c:8: RuntimeWarning: invalid value encountered in double_scalars
plt.plot(time_30s, cos_angles)
[<matplotlib.lines.Line2D at 0x10b3ae550>]