[R] Permutations
Rolf Turner
rolf at math.unb.ca
Wed Jul 14 15:06:41 CEST 2004
In respect of generating random ``restricted'' permutations, it
occurred to me as I was driving home last night .... If one is going
to invoke some kind of ``try again if it doesn't work procedure''
then one might as well keep it simple: Essentially use the rejection
method. Just generate a random permutation, and then check whether
it meets the restriction criterion. If yes, return that
permutation, if not, throw it away and try again.
This will definitely (???) genererate a genuinely random restricted
permutation. I figured that since a very large fraction of permutations
are acutally restricted permutions one wouldn't reject much of the
time, and so the rejection method should be pretty efficient.
I wrote the following code to do the work:
restr.perm2 <- function () {
#
okay <- function (x) {
m1 <- cbind(1:12,rep(1:4,rep(3,4)))
m2 <- m1[x,]
all((m2[,1] == m1[,1]) | (m2[,2] != m1[,2]))
}
#
repeat{
x <- sample(12,12)
if(okay(x)) return(x)
}
}
I'm bothered again, but: I did a little test to see how many tries
it would take on average. On 1000 attempts I got a mean of 8.44
tries, with a standard deviation of 7.7610 (standard error of the
mean = 0.2454). The value of 7.76 is roughly consistent with
sqrt(1-p.hat)/p.hat = 7.92 that one gets from the geometric
distribution.
This would indicate that the fraction of ``restricted'' permutations
is something like p.hat = 1/8.44 = 0.1184834. Which is a lot less
than the rough estimate of (4.7 x 10^6)/12! approx. = 0.9853 from
Robert Baskin's back-of-the-envelope calculations.
What's going on/wrong?
cheers,
Rolf Turner
rolf at math.unb.ca
More information about the R-help
mailing list