Confusion about NIPY hrf

nipy_hrf.py

import numpy as np
from scipy.stats import gamma

def gamma_params(peak, fwhm):
    """Return parameters to scipy.stats.gamma corresponding for an HRF shape.

    This was mostly copied from nipy.

    Parameters
    ----------
    peak : float
        time of response peak
    fwhm : float
        fwhm of response curve

    Returns
    -------
    shape : float
        shape parameter for gamma dist
    scale : float
        scale parameter for gamma dist

    """
    shape = np.power(peak / fwhm, 2) * 8 * np.log(2.0)
    scale = np.power(fwhm, 2) / peak / 8 / np.log(2.0)
    return shape, scale


if __name__ == "__main__":

    peak, fwhm = 5.4, 5.2 # Glover params
    shape, scale = gamma_params(peak, fwhm)
    gamma_rv = gamma(shape, scale=scale)
    x = np.linspace(0, 30, 1001)
    gamma_pdf = gamma_rv.pdf(x)
    print "Peak argument:", peak
    print "PDF peak: %.1f" % x[np.argmax(gamma_pdf)]