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

Ralf Gommers ralf.gommers@gmail....
Fri Mar 29 15:34:49 CDT 2013


On Fri, Mar 29, 2013 at 3:58 AM, <josef.pktd@gmail.com> wrote:

> 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.
>

Are you sure that works? Last time I checked we couldn't add docstrings to
f2py-wrapped functions with add_newdoc.

Ralf
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20130329/02dbc270/attachment.html 


More information about the SciPy-User mailing list