123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115 |
- function x = spm_invTcdf(F,v)
- % Inverse Cumulative Distribution Function (CDF) of Students t distribution
- % FORMAT x = spm_invTcdf(F,v)
- %
- % F - CDF (lower tail p-value)
- % v - degrees of freedom (v>0, non-integer d.f. accepted)
- % x - T-variate (Student's t has range (-Inf,Inf))
- %__________________________________________________________________________
- %
- % spm_invTcdf implements the inverse Cumulative Distribution of Students
- % t-distributions.
- %
- % Definition:
- %--------------------------------------------------------------------------
- % The Student's t-distribution with v degrees of freedom is defined for
- % positive integer v and real x. The Cumulative Distribution
- % Function (CDF) F(x) is the probability that a realisation of a
- % t-distributed random variable X has value less than x. F(x)=Pr{X<x}:
- % (See Evans et al., Ch37)
- %
- % This implementation is not restricted to whole (positive integer) df
- % v, rather it will compute for any df v>0.
- %
- % Variate relationships: (Evans et al., Ch37 & 7)
- %--------------------------------------------------------------------------
- % The Student's t distribution with 1 degree of freedom is the Standard
- % Cauchy distribution, which has a simple closed form CDF.
- %
- % For X a t-variate with v degrees of freedom, v/(v+X^2) has
- % distribution related to a Beta random variable with shape parameters
- % w/2 & 1/2, as described below.
- %
- % Algorithm:
- %--------------------------------------------------------------------------
- % Using the routine spm_invBcdf for the Beta distribution, with
- % appropriate parameters: The CDF of the Student's t-distribution with
- % v degrees of freedom is related to the incomplete beta function by:
- % Pr(|X|<x) = betainc(v/(v+x^2),v/2,1/2)
- % so
- % { betainc(v/(v+x^2),v/2,1/2) / 2 for x<0
- % F(x) = | 0.5 for x=0
- % { 1 - betainc(v/(v+x^2),v/2,1/2) / 2 for x>0
- %
- % See Abramowitz & Stegun, 26.5.27 & 26.7.1; Press et al., Sec6.4 for
- % definitions of the incomplete beta function. The relationship is
- % easily verified by substituting for v/(v+x^2) in the integral of the
- % incomplete beta function.
- %
- % References:
- %--------------------------------------------------------------------------
- % Evans M, Hastings N, Peacock B (1993)
- % "Statistical Distributions"
- % 2nd Ed. Wiley, New York
- %
- % Abramowitz M, Stegun IA, (1964)
- % "Handbook of Mathematical Functions"
- % US Government Printing Office
- %
- % Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
- % "Numerical Recipes in C"
- % Cambridge
- %
- %__________________________________________________________________________
- % Copyright (C) 1993-2011 Wellcome Trust Centre for Neuroimaging
- % Andrew Holmes
- % $Id: spm_invTcdf.m 4182 2011-02-01 12:29:09Z guillaume $
- %-Format arguments, note & check sizes
- %--------------------------------------------------------------------------
- if nargin<2, error('Insufficient arguments'), end
- ad = [ndims(F);ndims(v)];
- rd = max(ad);
- as = [ [size(F),ones(1,rd-ad(1))];...
- [size(v),ones(1,rd-ad(2))] ];
- rs = max(as);
- xa = prod(as,2)>1;
- if all(xa) && any(diff(as(xa,:)))
- error('non-scalar args must match in size');
- end
- %-Computation
- %--------------------------------------------------------------------------
- %-Initialise result to zeros
- x = zeros(rs);
- %-Only defined for F in [0,1] & strictly positive v.
- % Return NaN if undefined.
- md = ( F>=0 & F<=1 & v>0 );
- if any(~md(:))
- x(~md) = NaN;
- warning('Returning NaN for out of range arguments');
- end
- %-Special case: x is 0 when F=0.5, -Inf when F=0, +Inf when F=1
- x(md & F==0) = -Inf;
- x(md & F==1) = +Inf;
- %-Special case: Standard Cauchy distribution when v=1
- ml = ( md & v==1 ); if xa(1), mlF=ml; else mlF=1; end
- x(ml) = tan(pi*(F(mlF)-0.5));
- %-Compute where defined & not special cases
- Q = find( md & F>0 & F~=0.5 & F<1 & v~=1 );
- if isempty(Q), return, end
- if xa(1), QF=Q; else QF=1; end
- if xa(2), Qv=Q; else Qv=1; end
- %-Compute
- xQPos = F(QF)>0.5;
- bQ = spm_invBcdf(2*(xQPos -(xQPos*2-1).*F(QF)),v(Qv)/2,1/2);
- x(Q) = (xQPos*2-1) .* sqrt(v(Qv)./bQ -v(Qv));
|