/usr/share/pyshared/mvpa/measures/glm.py is in python-mvpa 0.4.8-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 | # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
#
# See COPYING file distributed along with the PyMVPA package for the
# copyright and license terms.
#
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
"""The general linear model (GLM)."""
__docformat__ = 'restructuredtext'
import numpy as N
from mvpa.measures.base import FeaturewiseDatasetMeasure
from mvpa.misc.state import StateVariable
class GLM(FeaturewiseDatasetMeasure):
"""General linear model (GLM).
Regressors can be defined in a design matrix and a linear fit of the data
is computed univariately (i.e. indepently for each feature). This measure
can report 'raw' parameter estimates (i.e. beta weights) of the linear
model, as well as standardized parameters (z-stat) using an ordinary
least squares (aka fixed-effects) approach to estimate the parameter
estimate.
The measure is reported in a (nfeatures x nregressors)-shaped array.
"""
pe = StateVariable(enabled=False,
doc="Parameter estimates (nfeatures x nparameters).")
zstat = StateVariable(enabled=False,
doc="Standardized parameter estimates (nfeatures x nparameters).")
def __init__(self, design, voi='pe', **kwargs):
"""
:Parameters:
design: array(nsamples x nregressors)
GLM design matrix.
voi: 'pe' | 'zstat'
Variable of interest that should be reported as feature-wise
measure. 'beta' are the parameter estimates and 'zstat' returns
standardized parameter estimates.
"""
FeaturewiseDatasetMeasure.__init__(self, **kwargs)
# store the design matrix as a such (no copying if already array)
self._design = N.asmatrix(design)
# what should be computed ('variable of interest')
if not voi in ['pe', 'zstat']:
raise ValueError, \
"Unknown variable of interest '%s'" % str(voi)
self._voi = voi
# will store the precomputed Moore-Penrose pseudo-inverse of the
# design matrix (lazy calculation)
self._inv_design = None
# also store the inverse of the inner product for beta variance
# estimation
self._inv_ip = None
def _call(self, dataset):
# just for the beauty of it
X = self._design
# precompute transformation is not yet done
if self._inv_design is None:
self._inv_ip = (X.T * X).I
self._inv_design = self._inv_ip * X.T
# get parameter estimations for all features at once
# (betas x features)
betas = self._inv_design * dataset.samples
# charge state
self.states.pe = pe = betas.T.A
# if betas and no z-stats are desired return them right away
if self._voi == 'pe' and not self.states.isEnabled('zstat'):
# return as (feature x beta)
return pe
# compute residuals
residuals = X * betas
residuals -= dataset.samples
# estimates of the parameter variance and compute zstats
# assumption of mean(E) == 0 and equal variance
# XXX next lines ignore off-diagonal elements and hence covariance
# between regressors. The humble being writing these lines asks the
# god of statistics for forgives, because it knows not what it does
diag_ip = N.diag(self._inv_ip)
# (features x betas)
beta_vars = N.array([ r.var() * diag_ip for r in residuals.T ])
# (parameter x feature)
zstat = pe / N.sqrt(beta_vars)
# charge state
self.states.zstat = zstat
if self._voi == 'pe':
# return as (feature x beta)
return pe
elif self._voi == 'zstat':
# return as (feature x zstat)
return zstat
# we shall never get to this point
raise ValueError, \
"Unknown variable of interest '%s'" % str(self._voi)
|