1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
|
#!/usr/bin/env python3
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
"""
In this example, we create a regression model for an event-related design in
which the response to an event at time T[i] is modeled as depending on the
amount of time since the last stimulus T[i-1]
"""
import numpy as np
import sympy
from nipy.algorithms.statistics.api import Formula, make_recarray
from nipy.modalities.fmri import hrf, utils
# Inter-stimulus intervals (time between events)
dt = np.random.uniform(low=0, high=2.5, size=(50,))
# Onset times from the ISIs
t = np.cumsum(dt)
# We're going to model the amplitudes ('a') by dt (the time between events)
a = sympy.Symbol('a')
linear = utils.define('linear', utils.events(t, dt, f=hrf.glover))
quadratic = utils.define('quad', utils.events(t, dt, f=hrf.glover, g=a**2))
cubic = utils.define('cubic', utils.events(t, dt, f=hrf.glover, g=a**3))
f1 = Formula([linear, quadratic, cubic])
# Evaluate this time-based formula at specific times to make the design matrix
tval = make_recarray(np.linspace(0,100, 1001), 't')
X1 = f1.design(tval, return_float=True)
# Now we make a model where the relationship of time between events and signal
# is an exponential with a time constant tau
l = sympy.Symbol('l')
exponential = utils.events(t, dt, f=hrf.glover, g=sympy.exp(-l*a))
f3 = Formula([exponential])
# Make a design matrix by passing in time and required parameters
params = make_recarray([(4.5, 3.5)], ('l', '_b0'))
X3 = f3.design(tval, params, return_float=True)
# the columns or d/d_b0 and d/dl
tt = tval.view(np.float64)
v1 = np.sum([hrf.glovert(tt - s)*np.exp(-4.5*a) for s,a in zip(t, dt)], 0)
v2 = np.sum([-3.5*a*hrf.glovert(tt - s)*np.exp(-4.5*a) for s,a in zip(t, dt)], 0)
V = np.array([v1,v2]).T
W = V - np.dot(X3, np.dot(np.linalg.pinv(X3), V))
np.testing.assert_almost_equal((W**2).sum() / (V**2).sum(), 0)
|
