Package mvpa :: Package misc :: Module fx

[hide private]
[frames] | no frames]

# Source Code for Module mvpa.misc.fx

      1  # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
      2  # vi: set ft=python sts=4 ts=4 sw=4 et:
      3  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
      4  #
      5  #   See COPYING file distributed along with the PyMVPA package for the
      6  #   copyright and license terms.
      7  #
      8  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
      9  """Misc. functions (in the mathematical sense)"""
     10
     11  __docformat__ = 'restructuredtext'
     12
     13  import numpy as N
     14
     15
     16 -def singleGammaHRF(t, A=5.4, W=5.2, K=1.0):
     17      """Hemodynamic response function model.
     18
     19      The version consists of a single gamma function (also see
     20      doubleGammaHRF()).
     21
     22      :Parameters:
     23        t: float
     24          Time.
     25        A: float
     26          Time to peak.
     27        W: float
     28          Full-width at half-maximum.
     29        K: float
     30          Scaling factor.
     31      """
     32      A = float(A)
     33      W = float(W)
     34      K = float(K)
     35      return K * (t / A) ** ((A ** 2) / (W ** 2) * 8.0 * N.log(2.0)) \
     36             * N.e ** ((t - A) / -((W ** 2) / A / 8.0 / N.log(2.0)))
     37
     38
     39 -def doubleGammaHRF(t, A1=5.4, W1=5.2, K1=1.0, A2=10.8, W2=7.35, K2=0.35):
     40      """Hemodynamic response function model.
     41
     42      The version is using two gamma functions (also see singleGammaHRF()).
     43
     44      :Parameters:
     45        t: float
     46          Time.
     47        A: float
     48          Time to peak.
     49        W: float
     50          Full-width at half-maximum.
     51        K: float
     52          Scaling factor.
     53
     54      Parameters A, W and K exists individually for each of both gamma
     55      functions.
     56      """
     57      A1 = float(A1)
     58      W1 = float(W1)
     59      K1 = float(K1)
     60      A2 = float(A2)
     61      W2 = float(W2)
     62      K2 = float(K2)
     63      return singleGammaHRF(t, A1, W1, K1) - singleGammaHRF(t, A2, W2, K2)
     64
     65
     66 -def leastSqFit(fx, params, y, x=None, **kwargs):
     67      """Simple convenience wrapper around SciPy's optimize.leastsq.
     68
     69      The advantage of using this wrapper instead of optimize.leastsq directly
     70      is, that it automatically constructs an appropriate error function and
     71      easily deals with 2d data arrays, i.e. each column with multiple values for
     72      the same function argument (`x`-value).
     73
     74      :Parameters:
     75        fx: functor
     76          Function to be fitted to the data. It has to take a vector with
     77          function arguments (`x`-values) as the first argument, followed by
     78          an arbitrary number of (to be fitted) parameters.
     79        params: sequence
     80          Sequence of start values for all to be fitted parameters. During
     81          fitting all parameters in this sequences are passed to the function
     82          in the order in which they appear in this sequence.
     83        y: 1d or 2d array
     84          The data the function is fitted to. In the case of a 2d array, each
     85          column in the array is considered to be multiple observations or
     86          measurements of function values for the same `x`-value.
     87        x: Corresponding function arguments (`x`-values) for each datapoint, i.e.
     88          element in `y` or columns in `y', in the case of `y` being a 2d array.
     89          If `x` is not provided it will be generated by `N.arange(m)`, where
     90          `m` is either the length of `y` or the number of columns in `y`, if
     91          `y` is a 2d array.
     92        **kwargs:
     93          All additonal keyword arguments are passed to `fx`.
     94
     95      :Returns:
     96        tuple: as returned by scipy.optimize.leastsq
     97          i.e. 2-tuple with list of final (fitted) parameters of `fx` and an
     98          integer value indicating success or failure of the fitting procedure
     99          (see leastsq docs for more information).
    100      """
    101      # import here to not let the whole module depend on SciPy
    102      from scipy.optimize import leastsq
    103
    104      y = N.asanyarray(y)
    105
    106      if len(y.shape) > 1:
    107          nsamp, ylen = y.shape
    108      else:
    109          nsamp, ylen = (1, len(y))
    110
    111      # contruct matching x-values if necessary
    112      if x is None:
    113          x = N.arange(ylen)
    114
    115      # transform x and y into 1d arrays
    116      if nsamp > 1:
    117          x = N.array([x] * nsamp).ravel()
    118          y = y.ravel()
    119
    120      # define error function
    121      def efx(p):
    122          err = y - fx(x, *p, **kwargs)
    123          return err
    124
    125      # do fit
    126      return leastsq(efx, params)
    127