[R] qt with ncp>37.62

[Charles C. Berry]

On Sat, 14 Jun 2008, data at tulipany.cz wrote:

> help(qt) states that:
> "ncp 	non-centrality parameter delta; currently except for rt(), only for
> abs(ncp) <= 37.62"
>
> so I would expect that calling qt with non-centrality parameter exceeding
> 37.62 should fail, instead e.g. calling
>
>> mapply(function(x) qt(p = 0.9, df = 55, ncp = x),35:45)
>
> gives:
>
> [1] 40.21448 41.35293 42.49164 43.68862 44.82945 45.97048 47.11170 48.25310
> [9] 49.39467 50.53639 51.67826
> Warning messages:
> 1: In qt(p = 0.9, df = 55, ncp = x) :
>  full precision was not achieved in 'pnt'
> 2: In qt(p = 0.9, df = 55, ncp = x) :
>  full precision was not achieved in 'pnt'
> 3: In qt(p = 0.9, df = 55, ncp = x) :
>  full precision was not achieved in 'pnt'
>
> so it seems calculation for (according to what is written in documentation)
> allowed values of ncp, i.e. in this case 35,36 and 37 is done and precision is
> checked, whereas calculation for the rest may be completely incorrect?
> Or was there any update of code (in pnt.c ?), allowing calculation of pt with
> higher ncp, not followed by documentation update?


The ultimate reference is the source code (see src/nmath/pnt.c), but
the lines below suggest that values smaller than 37.62 are not
guaranteed:


> qt( 0.9, df=55, ncp= 37 ) # WARNING ISSUED
[1] 42.49164
Warning message:
In qt(0.9, df = 55, ncp = 37) : full precision was not achieved in 'pnt'
>
>
> qt( 0.9, df=55, ncp= 38 ) # NO WARNING
[1] 43.68862
>

And the value for ncp=38 is too big, it seems:

table( rnorm(1000000,38)/sqrt(rchisq(1000000,55)/55)>qt( 0.9, df=55, ncp=38 ) )

and this is more pronounced at p == 0.9999.

So it appears that no checking is tried above abs(ncp)==37.62.

Further, there are discontinuities at abs(37.62) shown here:

library(splines)

plot( residuals( lm( qt(p = 0.5, df = 55, ncp=(-38):38) ~
 	bs(I((-38):38),knots=0,deg=1) ) ) )

and the value for ncp=37 shown here is accurate up to simulation error:

table( rnorm(1000000,37)/sqrt(rchisq(1000000,55)/55)>qt( 0.9, df=55, ncp=37 ) )

and seems to be about as accurate even at p == 0.9999

---

So, I guess what you get above 37.62 is based on some rougher
approximation and that it is not checked for accuracy.

HTH,

Chuck


>
> Thank you,
> Nikola Kaspříková
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

Charles C. Berry                            (858) 534-2098
                                             Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901