[SciPy-User] special's hidden treasures: rootfinding for distribution cdfs

josef.pktd@gmai... josef.pktd@gmai...
Thu Mar 28 21:58:14 CDT 2013


The fortran code for the distribution functions also has rootfinding build in.
We can solve for any value given the other values.

For example, solve for the non-centrality parameter for the
non-central t-distribution.
I never thought I would ever need that.

Is this the right function for special.nctdtrinc and friends? I didn't
manage to follow the trail.

\scipy\special\cdflib\cdftnc.f

C      SUBROUTINE CDFTNC( WHICH, P, Q, T, DF, PNONC, STATUS, BOUND )
C               Cumulative Distribution Function
C                  Non-Central T distribution
C
C                               Function
C
C     Calculates any one parameter of the noncentral t distribution give
C     values for the others.
C
C                               Arguments
C
C     WHICH --> Integer indicating which  argument
C               values is to be calculated from the others.
C               Legal range: 1..3
C               iwhich = 1 : Calculate P and Q from T,DF,PNONC
C               iwhich = 2 : Calculate T from P,Q,DF,PNONC
C               iwhich = 3 : Calculate DF from P,Q,T
C               iwhich = 4 : Calculate PNONC from P,Q,DF,T


I need iwhich=4

solving for the non-centrality parameter
>>> t = stats.nct._ppf(0.1, 9, 10)
>>> special.nctdtrinc(9, 0.1, t)
10.000000055648481

confidence interval for non-centrality parameter (AFAIU)

>>> special.nctdtrinc(24, [0.975, 0.025], 25)
array([ 17.68259,  32.26892])

upper value differs slightly from an R package (which has it's own
implementation)

cross check
>>> nc = special.nctdtrinc(24, 0.025, 25)
>>> stats.nct.cdf(25, 24, nc)
0.024999999999859579

a paper refers to a SAS function which seems to do the same
''''
following SAS syntax:
LowNC_CV = TNONCT(tobs., ν, 1 - α / 2),
'''

Background: If I understand this part of some papers correctly, then I can get
confidence intervals for effect sizes for statistical tests that are
based on the
t-distributions with non-central t-distribution under the alternative, with
essentially just a call to scipy.special.

------------------
It's nice to find the answer in scipy special when we need to look for
something we never knew that it existed.



**Scipy special needs your help.**
==========================

\scipy\special\add_newdocs.py      (maybe my checkout is a bit outdated)

<...>
add_newdoc("scipy.special", "ncfdtridfd",
    """
    """)

add_newdoc("scipy.special", "ncfdtridfn",
    """
    """)

add_newdoc("scipy.special", "ncfdtrinc",
    """
    """)

add_newdoc("scipy.special", "nctdtr",
    """
    """)

add_newdoc("scipy.special", "nctdtridf",
    """
    """)

add_newdoc("scipy.special", "nctdtrinc",
    """
    """)

add_newdoc("scipy.special", "nctdtrit",
    """
    """)
<...>

The information is available in the fortran docstrings. A pull request
that fills in some blanks is an easy way to become a scipy
contributer.


Josef


More information about the SciPy-User mailing list