algorithms.statistics.rft

Module: algorithms.statistics.rft

Inheritance diagram for nipy.algorithms.statistics.rft:

digraph inheritance17b0564be0 { rankdir=LR; size="8.0, 12.0"; "numpy.poly1d" [fontname="Vera Sans, DejaVu Sans, Liberation Sans, Arial, Helvetica, sans",fontsize=10,height=0.25,shape=box,style="setlinewidth(0.5)",tooltip="A one-dimensional polynomial class."]; "statistics.rft.ChiBarSquared" [URL="#nipy.algorithms.statistics.rft.ChiBarSquared",fontname="Vera Sans, DejaVu Sans, Liberation Sans, Arial, Helvetica, sans",fontsize=10,height=0.25,shape=box,style="setlinewidth(0.5)",target="_top"]; "statistics.rft.ChiSquared" -> "statistics.rft.ChiBarSquared" [arrowsize=0.5,style="setlinewidth(0.5)"]; "statistics.rft.ChiSquared" [URL="#nipy.algorithms.statistics.rft.ChiSquared",fontname="Vera Sans, DejaVu Sans, Liberation Sans, Arial, Helvetica, sans",fontsize=10,height=0.25,shape=box,style="setlinewidth(0.5)",target="_top",tooltip="EC densities for a Chi-Squared(n) random field."]; "statistics.rft.ECcone" -> "statistics.rft.ChiSquared" [arrowsize=0.5,style="setlinewidth(0.5)"]; "statistics.rft.ECcone" [URL="#nipy.algorithms.statistics.rft.ECcone",fontname="Vera Sans, DejaVu Sans, Liberation Sans, Arial, Helvetica, sans",fontsize=10,height=0.25,shape=box,style="setlinewidth(0.5)",target="_top",tooltip="EC approximation to supremum distribution of var==1 Gaussian process"]; "statistics.rft.IntrinsicVolumes" -> "statistics.rft.ECcone" [arrowsize=0.5,style="setlinewidth(0.5)"]; "statistics.rft.ECquasi" [URL="#nipy.algorithms.statistics.rft.ECquasi",fontname="Vera Sans, DejaVu Sans, Liberation Sans, Arial, Helvetica, sans",fontsize=10,height=0.25,shape=box,style="setlinewidth(0.5)",target="_top",tooltip="Polynomials with premultiplier"]; "numpy.poly1d" -> "statistics.rft.ECquasi" [arrowsize=0.5,style="setlinewidth(0.5)"]; "statistics.rft.FStat" [URL="#nipy.algorithms.statistics.rft.FStat",fontname="Vera Sans, DejaVu Sans, Liberation Sans, Arial, Helvetica, sans",fontsize=10,height=0.25,shape=box,style="setlinewidth(0.5)",target="_top",tooltip="EC densities for a F random field."]; "statistics.rft.ECcone" -> "statistics.rft.FStat" [arrowsize=0.5,style="setlinewidth(0.5)"]; "statistics.rft.Hotelling" [URL="#nipy.algorithms.statistics.rft.Hotelling",fontname="Vera Sans, DejaVu Sans, Liberation Sans, Arial, Helvetica, sans",fontsize=10,height=0.25,shape=box,style="setlinewidth(0.5)",target="_top",tooltip="Hotelling's T^2"]; "statistics.rft.ECcone" -> "statistics.rft.Hotelling" [arrowsize=0.5,style="setlinewidth(0.5)"]; "statistics.rft.IntrinsicVolumes" [URL="#nipy.algorithms.statistics.rft.IntrinsicVolumes",fontname="Vera Sans, DejaVu Sans, Liberation Sans, Arial, Helvetica, sans",fontsize=10,height=0.25,shape=box,style="setlinewidth(0.5)",target="_top",tooltip="Compute intrinsic volumes of products of sets"]; "statistics.rft.MultilinearForm" [URL="#nipy.algorithms.statistics.rft.MultilinearForm",fontname="Vera Sans, DejaVu Sans, Liberation Sans, Arial, Helvetica, sans",fontsize=10,height=0.25,shape=box,style="setlinewidth(0.5)",target="_top",tooltip="Maximize a multivariate Gaussian form"]; "statistics.rft.ECcone" -> "statistics.rft.MultilinearForm" [arrowsize=0.5,style="setlinewidth(0.5)"]; "statistics.rft.OneSidedF" [URL="#nipy.algorithms.statistics.rft.OneSidedF",fontname="Vera Sans, DejaVu Sans, Liberation Sans, Arial, Helvetica, sans",fontsize=10,height=0.25,shape=box,style="setlinewidth(0.5)",target="_top",tooltip="EC densities for one-sided F statistic"]; "statistics.rft.ECcone" -> "statistics.rft.OneSidedF" [arrowsize=0.5,style="setlinewidth(0.5)"]; "statistics.rft.Roy" [URL="#nipy.algorithms.statistics.rft.Roy",fontname="Vera Sans, DejaVu Sans, Liberation Sans, Arial, Helvetica, sans",fontsize=10,height=0.25,shape=box,style="setlinewidth(0.5)",target="_top",tooltip="Roy's maximum root"]; "statistics.rft.ECcone" -> "statistics.rft.Roy" [arrowsize=0.5,style="setlinewidth(0.5)"]; "statistics.rft.TStat" [URL="#nipy.algorithms.statistics.rft.TStat",fontname="Vera Sans, DejaVu Sans, Liberation Sans, Arial, Helvetica, sans",fontsize=10,height=0.25,shape=box,style="setlinewidth(0.5)",target="_top",tooltip="EC densities for a t random field."]; "statistics.rft.ECcone" -> "statistics.rft.TStat" [arrowsize=0.5,style="setlinewidth(0.5)"]; "statistics.rft.fnsum" [URL="#nipy.algorithms.statistics.rft.fnsum",fontname="Vera Sans, DejaVu Sans, Liberation Sans, Arial, Helvetica, sans",fontsize=10,height=0.25,shape=box,style="setlinewidth(0.5)",target="_top"]; }

Random field theory routines

The theoretical results for the EC densities appearing in this module were partially supported by NSF grant DMS-0405970.

Taylor, J.E. & Worsley, K.J. (2012). “Detecting sparse cone alternatives

for Gaussian random fields, with an application to fMRI”. arXiv:1207.3840 [math.ST] and Statistica Sinica 23 (2013): 1629-1656.

Taylor, J.E. & Worsley, K.J. (2008). “Random fields of multivariate

test statistics, with applications to shape analysis.” arXiv:0803.1708 [math.ST] and Annals of Statistics 36( 2008): 1-27

Classes

ChiBarSquared

class nipy.algorithms.statistics.rft.ChiBarSquared(dfn=1, search=[1])[source]

Bases: nipy.algorithms.statistics.rft.ChiSquared

__init__(dfn=1, search=[1])[source]

Initialize self. See help(type(self)) for accurate signature.

density(x, dim)

The EC density in dimension dim.

integ(m=None, k=None)
pvalue(x, search=None)
quasi(dim)

(Quasi-)polynomial parts of EC density in dimension dim

  • ignoring a factor of (2pi)^{-(dim+1)/2} in front.

ChiSquared

class nipy.algorithms.statistics.rft.ChiSquared(dfn, dfd=inf, search=[1])[source]

Bases: nipy.algorithms.statistics.rft.ECcone

EC densities for a Chi-Squared(n) random field.

__init__(dfn, dfd=inf, search=[1])[source]

Initialize self. See help(type(self)) for accurate signature.

density(x, dim)

The EC density in dimension dim.

integ(m=None, k=None)
pvalue(x, search=None)
quasi(dim)

(Quasi-)polynomial parts of EC density in dimension dim

  • ignoring a factor of (2pi)^{-(dim+1)/2} in front.

ECcone

class nipy.algorithms.statistics.rft.ECcone(mu=[1], dfd=inf, search=[1], product=[1])[source]

Bases: nipy.algorithms.statistics.rft.IntrinsicVolumes

EC approximation to supremum distribution of var==1 Gaussian process

A class that takes the intrinsic volumes of a set and gives the EC approximation to the supremum distribution of a unit variance Gaussian process with these intrinsic volumes. This is the basic building block of all of the EC densities.

If product is not None, then this product (an instance of IntrinsicVolumes) will effectively be prepended to the search region in any call, but it will also affect the (quasi-)polynomial part of the EC density. For instance, Hotelling’s T^2 random field has a sphere as product, as does Roy’s maximum root.

__init__(mu=[1], dfd=inf, search=[1], product=[1])[source]

Initialize self. See help(type(self)) for accurate signature.

pvalue(x, search=None)[source]
integ(m=None, k=None)[source]
density(x, dim)[source]

The EC density in dimension dim.

quasi(dim)[source]

(Quasi-)polynomial parts of EC density in dimension dim

  • ignoring a factor of (2pi)^{-(dim+1)/2} in front.

ECquasi

class nipy.algorithms.statistics.rft.ECquasi(c_or_r, r=0, exponent=None, m=None)[source]

Bases: numpy.poly1d

Polynomials with premultiplier

A subclass of poly1d consisting of polynomials with a premultiplier of the form:

(1 + x^2/m)^-exponent

where m is a non-negative float (possibly infinity, in which case the function is a polynomial) and exponent is a non-negative multiple of 1/2.

These arise often in the EC densities.

Examples

>>> import numpy
>>> from nipy.algorithms.statistics.rft import ECquasi
>>> x = numpy.linspace(0,1,101)
>>> a = ECquasi([3,4,5])
>>> a
ECquasi(array([3, 4, 5]), m=inf, exponent=0.000000)
>>> a(3) == 3*3**2 + 4*3 + 5
True
>>> b = ECquasi(a.coeffs, m=30, exponent=4)
>>> numpy.allclose(b(x), a(x) * numpy.power(1+x**2/30, -4))
True
__init__(c_or_r, r=0, exponent=None, m=None)[source]

Initialize self. See help(type(self)) for accurate signature.

denom_poly()[source]

Base of the premultiplier: (1+x^2/m).

Examples

>>> import numpy
>>> b = ECquasi([3,4,20], m=30, exponent=4)
>>> d = b.denom_poly()
>>> d
poly1d([ 0.03333333,  0.        ,  1.        ])
>>> numpy.allclose(d.c, [1./b.m,0,1])
True
change_exponent(_pow)[source]

Change exponent

Multiply top and bottom by an integer multiple of the self.denom_poly.

Examples

>>> import numpy
>>> b = ECquasi([3,4,20], m=30, exponent=4)
>>> x = numpy.linspace(0,1,101)
>>> c = b.change_exponent(3)
>>> c
ECquasi(array([  1.11111111e-04,   1.48148148e-04,   1.07407407e-02,
         1.33333333e-02,   3.66666667e-01,   4.00000000e-01,
         5.00000000e+00,   4.00000000e+00,   2.00000000e+01]), m=30.000000, exponent=7.000000)
>>> numpy.allclose(c(x), b(x))
True
compatible(other)[source]

Check compatibility of degrees of freedom

Check whether the degrees of freedom of two instances are equal so that they can be multiplied together.

Examples

>>> import numpy
>>> b = ECquasi([3,4,20], m=30, exponent=4)
>>> x = numpy.linspace(0,1,101)
>>> c = b.change_exponent(3)
>>> b.compatible(c)
True
>>> d = ECquasi([3,4,20])
>>> b.compatible(d)
False
>>>
deriv(m=1)[source]

Evaluate derivative of ECquasi

Parameters

m : int, optional

Examples

>>> a = ECquasi([3,4,5])
>>> a.deriv(m=2) 
ECquasi(array([6]), m=inf, exponent=0.000000)
>>> b = ECquasi([3,4,5], m=10, exponent=3)
>>> b.deriv()
ECquasi(array([-1.2, -2. ,  3. ,  4. ]), m=10.000000, exponent=4.000000)
property c

The polynomial coefficients

property coef

The polynomial coefficients

property coefficients

The polynomial coefficients

property coeffs

The polynomial coefficients

integ(m=1, k=0)

Return an antiderivative (indefinite integral) of this polynomial.

Refer to polyint for full documentation.

See also

polyint

equivalent function

property o

The order or degree of the polynomial

property order

The order or degree of the polynomial

property r

The roots of the polynomial, where self(x) == 0

property roots

The roots of the polynomial, where self(x) == 0

property variable

The name of the polynomial variable

FStat

class nipy.algorithms.statistics.rft.FStat(dfn, dfd=inf, search=[1])[source]

Bases: nipy.algorithms.statistics.rft.ECcone

EC densities for a F random field.

__init__(dfn, dfd=inf, search=[1])[source]

Initialize self. See help(type(self)) for accurate signature.

density(x, dim)

The EC density in dimension dim.

integ(m=None, k=None)
pvalue(x, search=None)
quasi(dim)

(Quasi-)polynomial parts of EC density in dimension dim

  • ignoring a factor of (2pi)^{-(dim+1)/2} in front.

Hotelling

class nipy.algorithms.statistics.rft.Hotelling(dfd=inf, k=1, search=[1])[source]

Bases: nipy.algorithms.statistics.rft.ECcone

Hotelling’s T^2

Maximize an F_{1,dfd}=T_dfd^2 statistic over a sphere of dimension k.

__init__(dfd=inf, k=1, search=[1])[source]

Initialize self. See help(type(self)) for accurate signature.

density(x, dim)

The EC density in dimension dim.

integ(m=None, k=None)
pvalue(x, search=None)
quasi(dim)

(Quasi-)polynomial parts of EC density in dimension dim

  • ignoring a factor of (2pi)^{-(dim+1)/2} in front.

IntrinsicVolumes

class nipy.algorithms.statistics.rft.IntrinsicVolumes(mu=[1])[source]

Bases: object

Compute intrinsic volumes of products of sets

A simple class that exists only to compute the intrinsic volumes of products of sets (that themselves have intrinsic volumes, of course).

__init__(mu=[1])[source]

Initialize self. See help(type(self)) for accurate signature.

MultilinearForm

class nipy.algorithms.statistics.rft.MultilinearForm(*dims, **keywords)[source]

Bases: nipy.algorithms.statistics.rft.ECcone

Maximize a multivariate Gaussian form

Maximized over spheres of dimension dims. See:

Kuri, S. & Takemura, A. (2001). ‘Tail probabilities of the maxima of multilinear forms and their applications.’ Ann. Statist. 29(2): 328-371.

__init__(*dims, **keywords)[source]

Initialize self. See help(type(self)) for accurate signature.

density(x, dim)

The EC density in dimension dim.

integ(m=None, k=None)
pvalue(x, search=None)
quasi(dim)

(Quasi-)polynomial parts of EC density in dimension dim

  • ignoring a factor of (2pi)^{-(dim+1)/2} in front.

OneSidedF

class nipy.algorithms.statistics.rft.OneSidedF(dfn, dfd=inf, search=[1])[source]

Bases: nipy.algorithms.statistics.rft.ECcone

EC densities for one-sided F statistic

See:

Worsley, K.J. & Taylor, J.E. (2005). ‘Detecting fMRI activation allowing for unknown latency of the hemodynamic response.’ Neuroimage, 29,649-654.

__init__(dfn, dfd=inf, search=[1])[source]

Initialize self. See help(type(self)) for accurate signature.

density(x, dim)

The EC density in dimension dim.

integ(m=None, k=None)
pvalue(x, search=None)
quasi(dim)

(Quasi-)polynomial parts of EC density in dimension dim

  • ignoring a factor of (2pi)^{-(dim+1)/2} in front.

Roy

class nipy.algorithms.statistics.rft.Roy(dfn=1, dfd=inf, k=1, search=[1])[source]

Bases: nipy.algorithms.statistics.rft.ECcone

Roy’s maximum root

Maximize an F_{dfd,dfn} statistic over a sphere of dimension k.

__init__(dfn=1, dfd=inf, k=1, search=[1])[source]

Initialize self. See help(type(self)) for accurate signature.

density(x, dim)

The EC density in dimension dim.

integ(m=None, k=None)
pvalue(x, search=None)
quasi(dim)

(Quasi-)polynomial parts of EC density in dimension dim

  • ignoring a factor of (2pi)^{-(dim+1)/2} in front.

TStat

class nipy.algorithms.statistics.rft.TStat(dfd=inf, search=[1])[source]

Bases: nipy.algorithms.statistics.rft.ECcone

EC densities for a t random field.

__init__(dfd=inf, search=[1])[source]

Initialize self. See help(type(self)) for accurate signature.

density(x, dim)

The EC density in dimension dim.

integ(m=None, k=None)
pvalue(x, search=None)
quasi(dim)

(Quasi-)polynomial parts of EC density in dimension dim

  • ignoring a factor of (2pi)^{-(dim+1)/2} in front.

fnsum

class nipy.algorithms.statistics.rft.fnsum(*items)[source]

Bases: object

__init__(*items)[source]

Initialize self. See help(type(self)) for accurate signature.

Functions

nipy.algorithms.statistics.rft.Q(dim, dfd=inf)[source]

Q polynomial

If dfd == inf (the default), then Q(dim) is the (dim-1)-st Hermite polynomial:

\[H_j(x) = (-1)^j * e^{x^2/2} * (d^j/dx^j e^{-x^2/2})\]

If dfd != inf, then it is the polynomial Q defined in [Worsley19947]

Parameters

dim : int

dimension of polynomial

dfd : scalar

Returns

q_poly : np.poly1d instance

References

Worsley19947(1,2)

Worsley, K.J. (1994). ‘Local maxima and the expected Euler characteristic of excursion sets of chi^2, F and t fields.’ Advances in Applied Probability, 26:13-42.

A ball-shaped search region of radius r.

nipy.algorithms.statistics.rft.binomial(n, k)[source]

Binomial coefficient

n!

c = ———

(n-k)! k!

Parameters

n : float

n of (n, k)

k : float

k of (n, k)

Returns

c : float

Examples

First 3 values of 4 th row of Pascal triangle

>>> [binomial(4, k) for k in range(3)]
[1.0, 4.0, 6.0]
nipy.algorithms.statistics.rft.mu_ball(n, j, r=1)[source]

j`th curvature of `n-dimensional ball radius r

Return mu_j(B_n(r)), the j-th Lipschitz Killing curvature of the ball of radius r in R^n.

nipy.algorithms.statistics.rft.mu_sphere(n, j, r=1)[source]

j`th curvature for `n dimensional sphere radius r

Return mu_j(S_r(R^n)), the j-th Lipschitz Killing curvature of the sphere of radius r in R^n.

From Chapter 6 of

Adler & Taylor, ‘Random Fields and Geometry’. 2006.

nipy.algorithms.statistics.rft.scale_space(region, interval, kappa=1.0)[source]

scale space intrinsic volumes of region x interval

See:

Siegmund, D.O and Worsley, K.J. (1995). ‘Testing for a signal with unknown location and scale in a stationary Gaussian random field.’ Annals of Statistics, 23:608-639.

and

Taylor, J.E. & Worsley, K.J. (2005). ‘Random fields of multivariate test statistics, with applications to shape analysis and fMRI.’

(available on http://www.math.mcgill.ca/keith

A spherical search region of radius r.

nipy.algorithms.statistics.rft.volume2ball(vol, d=3)[source]

Approximate volume with ball

Approximate intrinsic volumes of a set with a given volume by those of a ball with a given dimension and equal volume.