Here we try the simplest possible general linear model on some example data
%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
You are used to seeing these now:
import numpy as np
import matplotlib.pyplot as plt
import nibabel as nib
Let's get the data:
import os
HOME = os.path.expanduser('~')
DATA_PATH = os.path.join(HOME, 'data')
img_fname = os.path.join(DATA_PATH, 'ds105', 'sub001', 'BOLD', 'task001_run001', 'bold.nii.gz')
img = nib.load(img_fname)
n_vols = img.shape[-1]
n_vols
121
Now let's get the task. There were blocks of different types of objects in an on-off pattern. We'll pool the objects to get just 'on-off':
# from sub001/model/model001/onsets/task001_run001
block_starts_seconds = np.array([12, 48, 84, 120, 156, 192, 228, 264])
block_length_seconds = 22
TR = 2.5
Same in scans:
block_starts_scans = block_starts_seconds / TR
block_length_scans = block_length_seconds / TR
block_starts_scans
array([ 4.8, 19.2, 33.6, 48. , 62.4, 76.8, 91.2, 105.6])
Make a regressor. Let's give the HRF a delay of 2 scans or so.
delay = 2
regressor = np.zeros((n_vols,))
for start in block_starts_scans:
start_tr = round(start + delay)
end_tr = round(start + delay + block_length_scans)
regressor[start_tr:end_tr+1] = 1
print start_tr, end_tr
plt.plot(regressor)
7.0 16.0 21.0 30.0 36.0 44.0 50.0 59.0 64.0 73.0 79.0 88.0 93.0 102.0 108.0 116.0
[<matplotlib.lines.Line2D at 0x301bb90>]
Make the design:
X = np.column_stack((regressor, np.ones((n_vols,))))
Y = img.get_data()
slice_size = np.prod(Y.shape[:-1])
print slice_size
Y = np.reshape(Y, (slice_size, n_vols)).T
Y.shape
163840
(121, 163840)
X
array([[ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.]])
Fit the model!
B = np.linalg.pinv(X).dot(Y)
Pinv = np.linalg.pinv(X)
B = Pinv.dot(Y)
B = np.dot(Pinv, Y)
B.shape
(2, 163840)
Make B into an image again:
B = np.reshape(B, (2,) + img.shape[:-1])
B.shape
(2, 40, 64, 64)
Let's have a look at the image:
plt.imshow(B[0, :, :, 40], cmap='gray')
<matplotlib.image.AxesImage at 0x3140390>
B_selected = B[:, 19, 12, 40]
B_selected
array([ 14.78175313, 540.37209302])
Yselected = img.get_data()[19, 12, 40]
Yselected.shape
(121,)
Bselected = np.linalg.pinv(X).dot(Yselected)
Bselected
array([ 14.78175313, 540.37209302])
plt.plot(Yselected - np.mean(Yselected))
plt.plot(regressor * 40, 'r')
[<matplotlib.lines.Line2D at 0x3152c10>]
Xp = np.linalg.pinv(X)
Xp_again = np.linalg.inv(X.T.dot(X)).dot(X.T)
np.allclose(Xp - Xp_again, 0)
True
X2 = np.column_stack((X, X.dot([0.5, 0.5])))
# print X2
np.linalg.inv(X2.T.dot(X2))
--------------------------------------------------------------------------- LinAlgError Traceback (most recent call last) <ipython-input-34-56c3ce8ab61f> in <module>() ----> 1 np.linalg.inv(X2.T.dot(X2)) /home/jb/.local/lib/python2.7/site-packages/numpy/linalg/linalg.pyc in inv(a) 443 """ 444 a, wrap = _makearray(a) --> 445 return wrap(solve(a, identity(a.shape[0], dtype=a.dtype))) 446 447 /home/jb/.local/lib/python2.7/site-packages/numpy/linalg/linalg.pyc in solve(a, b) 326 results = lapack_routine(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0) 327 if results['info'] > 0: --> 328 raise LinAlgError('Singular matrix') 329 if one_eq: 330 return wrap(b.ravel().astype(result_t)) LinAlgError: Singular matrix
pinvX2 = np.linalg.pinv(X2)
print pinvX2.shape
(3, 121)