from __future__ import print_function
from __future__ import division
import numpy as np
import numpy.linalg as npl
import matplotlib.pyplot as plt
%matplotlib inline
Remember that the students-t statistic is:
$$ t = \frac{c^T \hat\beta}{\sqrt{\mathrm{var}(c^T \hat\beta)}} $$which works out to:
$$ t = \frac{c^T \hat\beta}{\sqrt{\hat{\sigma}^2 c^T (X^T X)^+ c}} $$where $\hat{\sigma}^2$ is our estimate of variance in the residuals, and $(X^T X)^+$ is the pseudo-inverse of $X^T X$.
from nipy.modalities.fmri.design import block_design, natural_spline, stack2designs
from nipy.modalities.fmri import hrf
times = np.arange(0, 240, 2) # TR onset times
T = len(times)
def random_onsets(times, p_thresh = 0.5):
is_onset = np.random.uniform(size=times.shape) > p_thresh
return times[is_onset]
isis = np.diff(random_onsets(times))
print(isis)
_ = plt.hist(isis)
[ 6 10 4 2 4 10 2 6 2 6 6 2 2 10 4 6 4 6 2 2 2 2 2 6 2 2 2 2 8 8 2 2 2 2 2 4 2 6 4 2 2 2 4 2 2 6 2 2 2 6 2 4 2 2 6 2 4 2 4 2 2 4 4 2]
onsets = random_onsets(times)
n_onsets = len(onsets)
offsets = onsets + 1.
block_spec = np.ones(n_onsets, dtype=([('start', float), ('end', float), ('type', int)]))
block_spec['start'] = onsets
block_spec['end'] = offsets
print(block_spec)
[(4.0, 5.0, 1) (6.0, 7.0, 1) (8.0, 9.0, 1) (10.0, 11.0, 1) (14.0, 15.0, 1) (16.0, 17.0, 1) (18.0, 19.0, 1) (20.0, 21.0, 1) (22.0, 23.0, 1) (24.0, 25.0, 1) (26.0, 27.0, 1) (28.0, 29.0, 1) (30.0, 31.0, 1) (32.0, 33.0, 1) (40.0, 41.0, 1) (48.0, 49.0, 1) (50.0, 51.0, 1) (60.0, 61.0, 1) (62.0, 63.0, 1) (64.0, 65.0, 1) (66.0, 67.0, 1) (68.0, 69.0, 1) (74.0, 75.0, 1) (76.0, 77.0, 1) (78.0, 79.0, 1) (82.0, 83.0, 1) (86.0, 87.0, 1) (90.0, 91.0, 1) (102.0, 103.0, 1) (106.0, 107.0, 1) (114.0, 115.0, 1) (118.0, 119.0, 1) (120.0, 121.0, 1) (122.0, 123.0, 1) (124.0, 125.0, 1) (126.0, 127.0, 1) (130.0, 131.0, 1) (132.0, 133.0, 1) (148.0, 149.0, 1) (158.0, 159.0, 1) (160.0, 161.0, 1) (166.0, 167.0, 1) (168.0, 169.0, 1) (172.0, 173.0, 1) (182.0, 183.0, 1) (184.0, 185.0, 1) (190.0, 191.0, 1) (192.0, 193.0, 1) (196.0, 197.0, 1) (204.0, 205.0, 1) (208.0, 209.0, 1) (210.0, 211.0, 1) (214.0, 215.0, 1) (216.0, 217.0, 1) (218.0, 219.0, 1) (222.0, 223.0, 1) (226.0, 227.0, 1) (230.0, 231.0, 1) (234.0, 235.0, 1) (236.0, 237.0, 1) (238.0, 239.0, 1)]
def random_design(times, duration=1):
onsets = random_onsets(times)
offsets = onsets + duration
n_onsets = len(onsets)
block_spec = np.ones(n_onsets, dtype=([('start', float), ('end', float), ('type', int)]))
block_spec['start'] = onsets
block_spec['end'] = offsets
return block_design(block_spec, times, hrfs=[hrf.spm])
X, contrasts = random_design(times)
plt.plot(times, X)
print(contrasts)
print(X.shape)
{'constant_0': array(1.0000000000000007)} (120,)
contrasts
{'constant_0': array(1.0000000000000007)}
spline = (natural_spline(times, order=3))
print(spline.shape)
(120, 4)
best_design = None
best_efficiency = 0
n_iters = 1000
efficiencies = np.zeros(n_iters)
splines = natural_spline(times)
for i in range(n_iters):
X_ev, contrasts_ev = random_design(times)
X, contrasts = stack2designs(X_ev, splines, contrasts_ev)
C = contrasts['constant_0'][:, None]
efficiency = 1. / C.T.dot(npl.pinv(X.T.dot(X)).dot(C))
efficiencies[i] = efficiency
if efficiency > best_efficiency:
best_efficiency = efficiency
best_design = X
stack2designs?
natural_spline?
import scipy.stats as sst
m = sst.scoreatpercentile(efficiencies, 95)
_ = plt.hist(efficiencies[efficiencies>m], bins=100)
efficiency
array([[ 109.45761821]])
best_efficiency
array([[ 227914.54247173]])
plt.plot(best_design[:,0])
[<matplotlib.lines.Line2D at 0x22cbccd0>]