Re: [R] qt with ncp>37.62

From: Charles C. Berry <cberry_at_tajo.ucsd.edu>
Date: Sun, 15 Jun 2008 11:37:14 -0700

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

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

Received on Sun 15 Jun 2008 - 18:42:07 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sun 15 Jun 2008 - 19:30:40 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive