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

josef.pktd@gmai... josef.pktd@gmai...
Fri Mar 29 15:55:31 CDT 2013


On Fri, Mar 29, 2013 at 4:34 PM, Ralf Gommers <ralf.gommers@gmail.com> wrote:
>
>
>
>
> 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.
>


Now I'm not sure, and I never tried.

I was just jumping to conclusions based on superficially following
https://github.com/scipy/scipy/pull/459

Josef



>
>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>


More information about the SciPy-User mailing list