Function signature

glm

hrf_estimation.glm(conditions, onsets, TR, Y, drifts=None, basis='3hrf', mode='r1glm', hrf_length=20, oversample=5, rtol=1e-08, verbose=False, maxiter=500, callback=None, method='L-BFGS-B', n_jobs=1, hrfs=None, return_design_matrix=False, bounds=True)

Perform a GLM from BOLD signal, given the conditons, onset, TR (repetition time of the scanner) and the BOLD signal.

This method is able to fir a variety of models, available through the mode keyword. These are:

• glm: standard GLM
• glms: GLM with separate designs
• r1glm: Rank-1 GLM
• r1glms: Rank-1 GLM with separate designs

basis:

• hrf: single element basis
• 3hrf: basis with 3 elements
• fir: basis with hrf_length elements (in multiples of TR)

Note the output parameters need are not normalized. Rank-1 models are specified up to a constant term between the betas and the HRF. This implies that some normalization must be done prior to interpreting the activation coefficients. Typically the HRF is normalized to have unit amplitude and to correlate positively with a reference HRF.

Parameters:

conditions: array-like, shape (n_trials) :

array of conditions

onsets: array-like, shape (n_trials) :

array of onsets

TR: float :

Repetition Time, the delay between two succesive aquisitions of the same image.

Y : array-like, shape (n_scans, n_voxels)

Time-series vector.

mode: {‘r1glm’, ‘r1glms’, ‘glms’, ‘glm’} :

Different GLM models.

rtol : float

Relative tolerance for stopping criterion.

maxiter : int

maximum number of iterations

verbose : {0, 1, 2}

Different levels of verbosity

n_jobs: int :

Number of CPUs to use. Use -1 to use all available CPUs.

method: {‘L-BFGS-B’, ‘TNC’} :

Different algorithmic solvers, only used for ‘r1*’ modes. All should yield the same result but their efficiency might vary.

Returns:

U : array

Estimated HRF. Will be of shape (basis_len, n_voxels) for rank-1 methods and of (basis_len, n_conditions, n_voxels) for the other methods.

V : array, shape (p, n_voxels)

Estimated activation coefficients (beta-map).

dmtx: array, :

Design matrix. Only returned if return_design_matrix=True

rank_one

hrf_estimation.rank_one(X, y, n_basis, w_i=None, drifts=None, callback=None, maxiter=500, method='L-BFGS-B', rtol=1e-06, verbose=0, mode='r1glm', hrfs=None, basis=None, bounds=True)

Estimates a R1-GLM model with a given design matrix.

This methods solves a problem of the form

argmin_{hrf, betas} ||y - X vec(hrf betas.T)||^2

for given y and X

Parameters:

X : array-like, shape (n_scans, n_regressors)

Design matrix. Note that the number of columns in the design matrix must be a multiple of the number of basis elements (n_basis)

y_i: array-like, shape (n_scans, n_voxels) :

BOLD signal.

n_basis : int

Number of basis elements in the HRF.

w_i : array-like

initial point.

method: {‘L-BFGS-B’, ‘TNC’} :

Optimization algorithm.

bounds: {True, False} :

If True, constraints the estimated HRF (only in the case of basis=‘2hrf’ or basis=‘3hrf’) so that the second and third derivative have absolute value < 1.

Returns:

U : array

Estimated HRFs

V : array

Estimated activation coefficients

spmt

hrf_estimation.hrf.spmt(t)

SPM canonical HRF, HRF values for time values t

This is the canonical HRF function as used in SPM. It has the following defaults:

defaults (seconds)

delay of response (relative to onset) 6 delay of undershoot (relative to onset) 16 dispersion of response 1 dispersion of undershoot 1 ratio of response to undershoot 6 onset (seconds) 0 length of kernel (seconds) 32

dspmt

hrf_estimation.hrf.dspmt(t)

SPM canonical HRF derivative, HRF derivative values for time values t

This is the canonical HRF derivative function as used in SPM.

It is the numerical difference of the HRF sampled at time t minus the values sampled at time t -1

ddspmt

hrf_estimation.hrf.ddspmt(t)

SPM canonical HRF dispersion derivative, values for time values t

This is the canonical HRF dispersion derivative function as used in SPM.

It is the numerical difference between the HRF sampled at time t, and values at t for another HRF shape with a small change in the peak dispersion parameter (peak_disp in func:spm_hrf_compat).

savgol_filter

hrf_estimation.savitzky_golay.savgol_filter(x, window_length, polyorder, deriv=0, delta=1.0, axis=-1, mode='interp', cval=0.0)

Apply a Savitzky-Golay filter to an array.

This is a 1-d filter. If x has dimension greater than 1, axis determines the axis along which the filter is applied.

Parameters:

x : array_like

The data to be filtered. If x is not a single or double precision floating point array, it will be converted to type numpy.float64 before filtering.

window_length : int

The length of the filter window (i.e. the number of coefficients). window_length must be a positive odd integer.

polyorder : int

The order of the polynomial used to fit the samples. polyorder must be less than window_length.

deriv : int, optional

The order of the derivative to compute. This must be a nonnegative integer. The default is 0, which means to filter the data without differentiating.

delta : float, optional

The spacing of the samples to which the filter will be applied. This is only used if deriv > 0. Default is 1.0.

axis : int, optional

The axis of the array x along which the filter is to be applied. Default is -1.

mode : str, optional

Must be ‘mirror’, ‘constant’, ‘nearest’ or ‘interp’. This determines the type of extension to use for the padded signal to which the filter is applied. When mode is ‘constant’, the padding value is given by cval. See the Notes for more details on ‘mirror’, ‘constant’ and ‘nearest’. When the ‘interp’ mode is selected (the default), no extension is used. Instead, a degree polyorder polynomial is fit to the last window_length values of the edges, and this polynomial is used to evaluate the last window_length // 2 output values.

cval : scalar, optional

Value to fill past the edges of the input if mode is ‘constant’. Default is 0.0.

Returns:

y : ndarray, same shape as x

The filtered data.

See also

savgol_coeffs

Notes

Details on the mode options:

‘mirror’:

Repeats the values at the edges in reverse order. The value closest to the edge is not included.

‘nearest’:

The extension contains the nearest input value.

‘constant’:

The extension contains the value given by the cval argument.

For example, if the input is [1, 2, 3, 4, 5, 6, 7, 8], and window_length is 7, the following shows the extended data for the various mode options (assuming cval is 0):

mode       |   Ext   |         Input          |   Ext
-----------+---------+------------------------+---------
'mirror'   | 4  3  2 | 1  2  3  4  5  6  7  8 | 7  6  5
'nearest'  | 1  1  1 | 1  2  3  4  5  6  7  8 | 8  8  8
'constant' | 0  0  0 | 1  2  3  4  5  6  7  8 | 0  0  0

Examples

>>> np.set_printoptions(precision=2)  # For compact display.
>>> x = np.array([2, 2, 5, 2, 1, 0, 1, 4, 9])

Filter with a window length of 5 and a degree 2 polynomial. Use the defaults for all other parameters.

>>> y = savgol_filter(x, 5, 2)
array([ 1.66,  3.17,  3.54,  2.86,  0.66,  0.17,  1.  ,  4.  ,  9.  ])

Note that the last five values in x are samples of a parabola, so when mode=’interp’ (the default) is used with polyorder=2, the last three values are unchanged. Compare that to, for example, mode=’nearest’:

>>> savgol_filter(x, 5, 2, mode='nearest')
array([ 1.74,  3.03,  3.54,  2.86,  0.66,  0.17,  1.  ,  4.6 ,  7.97])

create_design_matrix

hrf_estimation.utils.create_design_matrix(conditions, onsets, TR, n_scans, basis='3hrf', oversample=10, hrf_length=32)

Parameters:

conditions: list of conditions :

onset: list of onsets :

TR: float :

repetition time

n_scans: number of scans :

basis: one of {‘fir’, ‘3hrf’, ‘2hrf’, ‘hrf’} :

Returns:

design_matrix :

basis :