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