[R] family question

[Douglas Bates]

In the continuing saga of simulating family sizes from Troels Ring's
demographic problem, I offer a somewhat improved function.  In the previous
function if the process did not terminate in the current sample I
would append a new sample to it and recalculate the cumulative sum.
Obviously this is a waste because you already know that the process
will not terminate in the current sample and you know the total from
this sample.  The new function keeps track of the current total (tot)
and the family size to this point (offset) and only calculates
cumulative sums for the new part.

This function increases the size of chunks on which it is working up
to a threshold (2^20 by default) then fixes the chunk size and
continues until the total size exceeds another threshold (2^25 by
default).  I ran this simulation on cvs.r-project.org, a 400 MHz
Pentium III with 768 MB memory running Linux.  I could simulate
10000 samples in about 11 minutes.  The largest value is over
18,000,000 and there is one NA, which would have been over 32,000,000.

> fam
function(basesize = 16, maxvec = 2^20, maxoffset = 2^25)
{
    ## simulate the family sizes in Troels Ring's problem
    samp <- ifelse(runif(basesize) > 0.5, 1, -1)
    tot <- offset <- 0
    repeat {
        cs <- tot + cumsum(samp)
        moreboys <- cs > 0
        if (any(moreboys)) return(min((offset + seq(along = moreboys))[moreboys]))
        offset <- offset + length(samp)
        if (offset > maxoffset) return(NA)
        tot <- cs[length(cs)]
        samp <- ifelse(runif(min(2*length(samp), maxvec)) > 0.5, 1, -1)
    }
}
> famsamp
function(n = 10000)
{
   val <- integer(n)
   for (i in seq(along = val)) val[i] <- fam()
   val
}
> system.time(fs <- famsamp())
[1] 673.45   1.12 752.76   0.00   0.00
> summary(fs)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
       1        1        1     5091        9 18180000        1 
> sort(fs[fs > 100000])
 [1]   137031   156689   168471   190025   220981   221267   285271   379793
 [9]   423175   549397   571879   741895   766307   958423  1819621  1850911
[17]  2024143  2110585  2749441  2860301 11046833 18181117


-- 
Douglas Bates                            bates at stat.wisc.edu
Statistics Department                    608/262-2598
University of Wisconsin - Madison        http://www.stat.wisc.edu/~bates/
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._