Re: [R] Permutations

From: Rolf Turner <rolf_at_math.unb.ca>
Date: Wed 14 Jul 2004 - 23:06:41 EST


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@math.unb.ca

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Wed Jul 14 23:16:45 2004

This archive was generated by hypermail 2.1.8 : Fri 18 Mar 2005 - 02:36:05 EST