# Re: [R] Random numbers from noncentral t-distribution

From: Long Qu <longeregnol_at_yahoo.com.cn>
Date: Sat 01 Jul 2006 - 02:46:14 EST

Thank you very much for your kind reply. It solved the problem of rt( ). :D

But it seems that the qt( ) also have problems:   I modified the rt( ) function as you suggested,   rt <- function (n, df, ncp = 0)
{
if (ncp == 0)
.Internal(rt(n, df))
else rnorm(n, ncp)/sqrt(rchisq(n, df)/df) }

Then I increase the number of random variables to 10000, and made a QQ-plot:   > qqplot(rt(10000,df=20,ncp=3),qt(runif(10000),df=20,ncp=3))

I've got some spurious points at lower-left corner. It seems that qt( ) results were truncated.

I also tried this with another df and ncp:   > pt(-.75,df=2,ncp=1)
 0.05726429
> sum(qt(1:10000/10001,df=2,ncp=1)< -.75)
 0
where I'd expected the last number should be > 550 or so, not 0.

Thanks again, the modified rt( ) is now OK for my work.   Long

Thomas Lumley <tlumley@u.washington.edu> wroteㄩ   On Fri, 30 Jun 2006, Long Qu wrote:

> Hi there:
>
> I'd thought these two versions of noncentral t-distribution are essentially the same:
> > qqplot(rt(1000,df=20,ncp=3),qt(runif(1000),df=20,ncp=3))
>
> But, the scales of the x-axis and the y-axis are quite different according to the QQ-plot.
>
> Did I make any mistakes somewhere?
>

No, I think we did.

We have
> rt

function (n, df, ncp = 0)
{
if (ncp == 0)
.Internal(rt(n, df))
else rnorm(n, ncp)/(rchisq(n, df)/sqrt(df)) }

and the rchisq() in the denominator should be inside the sqrt().

-thomas

ツ蛁捇誥轎煤蚘眊-3.5G⺍講ㄛ20M蜇璃ㄐ

[[alternative HTML version deleted]]

R-help@stat.math.ethz.ch mailing list