[Rd] noncentral F-distributed random numbers (PR#9055)

From: <longor_at_iastate.edu>
Date: Sat 01 Jul 2006 - 19:13:47 GMT


Full_Name: Long Qu
Version: 2.3.1
OS: Windows XP
Submission from: (NULL) (64.113.93.235)

The QQ-plot of two versions of simulating noncentral F-distributed random numbers has quite different scales:
> qqplot(rf(1000,2,15,3),qf(runif(1000),2,15,3))

The rf() function reads:
> rf

function (n, df1, df2, ncp = 0)
{

    if (ncp == 0)

        .Internal(rf(n, df1, df2))
    else rchisq(n, df1, ncp = ncp)/rchisq(n, df2) }
<environment: namespace:stats>

where I believe both the numerator and the denominator should be divided by their corresponding degrees of freedom.

My suggested (slighly augmented) version is:
> rf=

function (n, df1, df2, ncp1 = 0, ncp2=0) {

    if (ncp1 == 0 && ncp2==0)

        .Internal(rf(n, df1, df2))
    else if (ncp2==0)

          (rchisq(n, df1, ncp = ncp1)/df1)/(rchisq(n, df2)/df2)     else (rchisq(n, df1, ncp = ncp1)/df1)/(rchisq(n, df2, ncp=ncp2)/df2)

}

which incorporated both singly and doubly noncentral F-distribution.

> version

               _                         
platform       i386-pc-mingw32           
arch           i386                      
os             mingw32                   
system         i386, mingw32             
status                                   
major          2                         
minor          3.1                       
year           2006                      
month          06                        
day            01                        
svn rev        38247                     
language       R                         
version.string Version 2.3.1 (2006-06-01)

Regards,
Long Qu



R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Sun Jul 02 05:16:15 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Mon 03 Jul 2006 - 20:24:51 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.