[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