probability


30
Aug 10

The Chosen One

Toss one hundred different balls into your basket. Shuffle them up and select one with equal probability amongst the balls. That ball you just selected, it’s special. Before you put it back, increase its weight by 1/100th. Then put it back, mix up the balls and pick again. If you do this enough, at some point there will be a consistent winner which begins to stand out.

The graph above shows the results of 1000 iterations with 20 balls (each victory increases the weight of the winner by 5%). The more balls you have, the longer it takes before a clear winner appears. Here’s the graph for 200 balls (0.5% weight boost for each victory).

As you can see, in this simulation it took about 85,000 iterations before a clear winner appeared.

I contend that as the number of iterations grows, the probability of seeing a Chosen One approaches unity, no matter how many balls you use. In other words, for any number of balls, a single one of them will eventually see its relative weight, compared to the others, diverge. Can you prove this is true?

BTW this is a good Monte Carlo simulation of the Matthew Effect (no relation).

Here is the code in R to replicate:

numbItems = 200
items = 1:numbItems
itemWeights = rep(1/numbItems,numbItems) # Start out uniform
iterations = 100000
itemHistory = rep(0,iterations)
 
for(i in 1:iterations) {
	chosen = sample(items, 1, prob=itemWeights)
	itemWeights[chosen] = itemWeights[chosen] + (itemWeights[chosen] * (1/numbItems))
	itemWeights = itemWeights / sum(itemWeights) # re-Normalze
	itemHistory[i] = chosen
}
 
plot(itemHistory, 1:iterations, pch=".", col="blue")

After many trials using a fixed large number of balls and iterations, I found that the moment of divergence was amazingly consistent. Do you get the same results?


8
Aug 10

Seeing angels in the architecture

Sorry for the long delay between posts; I was temporarily sucked in to the infinite. While doing some reading about set theory (foundational stuff for probability and, in fact, all of mathematics), I veered off into the infinite and had a hard time climbing back out. I’m guessing you already know most of the basics about sets: compliments and unions and intersections. You may even know some of the stranger parts, like G. Cantor’s cascading crescendo of cardinalities. But knowing those in a cursory way (and really, that’s all a work-a-day statistician or even probabilist needs) isn’t the same as really exploring them.

Looking up again now after several weeks, I feel like I’ve traveled three levels deep in a dream, lost in a purgatory I could only escape by answering questions likeĀ  “Is a line made up of points, or does it have points?”, “Is it possible to count what you cannot fully name”, and “In an unbounded universe, is the compliment of the compliment of an object the same as the original object?”. I know, I know. I should have taken that left back at Albuquerque, I shouldn’t have swallowed the red pill. Still, it’s been an interesting trip to say the least, and I feel like I may now be coming back up the the surface, a little bit wiser and a lot more confused than when I began.

Meanwhile, I’ve added a couple items to the “Manifesto” and, The Architect permitting, will be posting a theory on Types of Randomness soon. Post should take between 1 and 10 days to complete, with 95% confidence. Hum… better make that an 80% confidence interval, I still haven’t wrapped my head around the whole idea of forcing.


19
Jul 10

R: Clash of the cannon cycles

Imagine a unit square. Every side has length 1, perfectly square. Now imagine this square was really a fence, and you picked two spots at random along the fence, with uniform probability over the length of the fence. At each of these two locations, set down a special kind of cannon. Just like the light cycles from Tron, these cannons leave trails of color. To aim each cannon, pick another random point from one of the three other sides of the fence, and aim for that point.

Sometimes there will be a collision within the square, other times no. The image at top shows the results of five trials. The red dots are where the trails from a pair of cannons collided. My burning question: What is the distribution for these dots? Before reading on, try to make a guess. Where will collisions be most likely to happen?

Somewhere in the world, there lives a probabilist who could come up with a formula for the exact distribution in an hour, but that person doesn’t live in my house, so I took the Monte Carlo approach, coded in R:

# Functions to pick two points, not on the same side:
m2pt <- function(m) {
	if(m <1) {
		myPoint = c(5,m,0)
	} else if (m < 2) {
		myPoint = c(6,1,m %% 1) 
	} else if (m < 3) {
		myPoint = c(7,1-(m %% 1),1)
	} else {
		myPoint = c(8,0,1-(m %% 1))
	}
	return(myPoint)
}		
 
get2pts <- function() {
	pt1 = m2pt(runif(1,0,4))
	pt2 = m2pt(runif(1,0,4))
 
	# Make sure not both on the same sides. If so, keep trying
	while(pt1[1] == pt2[1]) {
		pt2 = m2pt(runif(1,0,4))
	}
	return(c(pt1[2:3],pt2[2:3]))
}
 
# Optional plot of every cannon fire line. Not a good idea for "iters" more than 100
#windows()
#plot(0,0,xlim=c(0,1),ylim=c(0,1),col="white")		
 
# How many times to run the experiment
iters = 10000
 
# Track where the intersections occur
interx = c()
intery = c()
 
for(i in 1:iters) {
	can1 = get2pts()
	can2 = get2pts()
 
	# Optional plot of every cannon fire line. Not a good idea for "iters" more than 100 
	#points(c(can1[1],can1[3]),c(can1[2],can1[4]),pch=20,col="yellow")
	#segments(can1[1],can1[2],can1[3],can1[4],pch=20,col="yellow",lwd=1.5)
	#points(c(can2[1],can2[3]),c(can2[2],can2[4]),pch=20,col="blue")
	#segments(can2[1],can2[2],can2[3],can2[4],pch=20,col="blue",lwd=1.5)
 
	# See if there is a point of intersection, find it. 
	toSolve = matrix(c( can1[3]-can1[1], can2[3]-can2[1], can1[4]-can1[2], can2[4]-can2[2]),byrow=T,ncol=2)
	paras = solve(toSolve, c( can2[1]-can1[1], can2[2]-can1[2]))
 
	solution = c(can1[1] + paras[1]*(can1[3]-can1[1]), can1[2] + paras[1]*(can1[4]-can1[2]))
 
	# Was the collision in the square
	if(min(solution) > 0 && max(solution) < 1) {
		# Optional plot of red dots
		# points(solution[1],solution[2],pch=20,col="red",cex=1.5)
 
		# if this intersection is in square, plot it, add it to list of intersections
		interx = c(interx, solution[1])
		intery = c(intery, solution[2])
	}
}
 
windows()
plot(interx, intery, pch=20,col="blue",xlim=c(0,1),ylim=c(0,1))

After carefully writing and debugging much more code than I expected, I ran a trial with several thousand cannon fires and plotted just the collisions. Here is what I saw:

Looks pretty uniform, doesn’t it? If it is, I will have gone a very long way just to replicate the bi-variate uniform distribution. My own original guess was that most collisions, if they happened in the square, would be towards the middle. Clearly this wasn’t the case. Looks can be deceiving, though, so I checked a histogram of the x’s (no need to check the y’s, by symmetry they have the same distribution):

Very interesting, no? The area near the edges appears more likely to have a collision, with an overall rounded bowl shape to the curve. The great thing about Monte Carlo simulations is that if something unexpected happens, you can always run it again with more trials. Here I changed “iters” to 100,000, ran the code again, and plotted the histogram.

hist(interx, breaks=100, col="blue",xlab="x",main="Histogram of the x's")

Now its clear that the distribution spikes way up near the edges, and appears to be essentially flat for most of the middle area. It seems like it may even go up slightly at the very middle. Just to be sure, I ran a trial with one million iterations:

Now it definitely looks like a small upward bulge in the middle, though to be sure I would have to do run some statistical tests or use an even larger Monte Carlo sample, and given how inefficient my code is, that could take the better part of a week to run. So for today I’ll leave it at that.

One final statistic of note: During my run of one million iterations, 47.22% of all collisions happened inside the box. What do you think, is the true, theoretical ratio of collisions within the box a rational number?


10
Jul 10

Photo without caption


9
Jul 10

100 Prisoners, 100 lines of code

In math and economics, there is a long, proud history of placing imaginary prisoners into nasty, complicated scenarios. We have, of course, the classic Prisoner’s Dilemma, as well as 100 prisoners and a light bulb. Add to that list the focus of this post, 100 prisoners and 100 boxes.

In this game, the warden places 100 numbers in 100 boxes, at random with equal probability that any number will be in any box. Each convict is assigned a number. One by one they enter the room with the boxes, and try to find their corresponding number. They can open up to 50 different boxes. Once they either find their number or fail, they move on to a different room and all of the boxes are returned to exactly how they were before the prisoner entered the room.

The prisoners can communicate with each other before the game begins, but as soon as it starts they have no way to signal to each other. The warden is requiring that all 100 prisoners find their numbers, otherwise she will force them to listen to hundreds of hours of non-stop, loud rock musician interviews. Can they avoid this fate?

The first thing you might notice is that if every prisoner opens 50 boxes at random, they will have a 0.5 probability of finding their number. The chances that all of them will find their number is (\frac{1}2)^{100}, which is approximately as rare as finding a friendly alien with small eyes. Can they do better?

Of course they can. Otherwise I wouldn’t be asking the question, right? Before I explain how, and go into a Monte Carlo simulation in R, you might want to think about how they can do it. No Googling!

All set? Did you find a better way? The trick should be clear from the code below, but if not skip on to the explanation.

# How many times should we run this experiment?
iters = 1000
results = rep(0,iters)
 
for(i in 1:iters) {
	# A random permutation:
	boxes = sample(1:100,100)
 
	# Labels for our prisoners
	prisoners = 1:100
 
	# Track how many "winners" we have
	foundIt = 0
 
	# Main loop over the prisoners
	for(prisoner in prisoners) {
 
		# Track the prisoners path
		path = c(prisoner)
 
		tries = 1
 
		# Look first in the box that matches your own number
		inBox = boxes[prisoner]
 
		while(tries &lt; 50) { 			 			path = c(path, inBox) 			 			if(inBox == prisoner) { 				#cat("Prisoner", prisoner, "found her number on try", tries, "\n") 				foundIt = foundIt + 1 				break; 			} else { 				# Follow that number to the next box 				inBox = boxes[inBox] 			} 			tries = tries+1 		} 		 		# cat("Prisoner", prisoner, "took path", paste(path, collapse=" -&gt; "), "\n")
	}
 
	# How many prisoners found their numbers?
	cat("A total of", foundIt, "prisoners found their numbers.\n")
	flush.console()
	results[i] = foundIt
}
 
hist(results, breaks=100, col="blue")

Here is what one of my plots looked like after running the code:

Out of the 1000 times I ran the experiment, on 307 occasions every single prisoner found his number. The theoretical success rate is about 31%. So, if it’s not clear from the code, what was the strategy employed by the prisoners and how does it work?

One way to look at the distribution of numbers in boxes is to see it as a permutation of the numbers from 1 to 100. Each permutation can be partitioned into what are called cycles. A cycle works like this: pick any number in your permutation. Let’s say it’s 23. Then you look at the number the 23rd place (ie the number in the 23rd box, counting from the left). If that number is 16, you look at the number in the 16th place. If that number is 87, go open box number 87 and follow that number. Eventually, the box you open up will have the number that brings you back to where you started, completing the cycle. Different permutations have different cycles.

The key for the prisoner is that by starting with the box that is the same place from the left as his number, and by following the numbers in the boxes, the prisoner guarantees that if he is in a cycle of length less than 50, he will eventually open the box with his number in it, which would complete the cycle he began. One way to envision cycles of different lengths is to think about the extreme cases. If a particular permutation shifted every single number over one to the left (and wrapped number 1 onto the end), you would have a single cycle of length 100. Box 1 would contain number 2, box 2 number 3 and so on. On the other hand, if a permutation flipped every pair of consecutive numbers, you would have 50 cycles, each of length 2: box 1 would have number 2, box 2 would have number 1. Of course if your permutation doesn’t change anything you have 100 cycles of length 1.

As you can see from the histogram, when using this strategy you can never have between 50 and 100 winning prisoners. Anytime you have a single cycle of length greater than 50, for example 55, then all 55 prisoners who start on that cycle will fail to find their number. If no cycles are longer than 50, everyone wins. Just how rare are different cycles of different lengths? For the math behind that check out this excellent explanation by Peter Taylor of Queen’s University.

Before moving on I wanted to visualize these cycles. Try running the code below:

# Unit circle
plot(0,0,xlim=c(-1,1),ylim=c(-1,1),col="white",ann=FALSE, xaxt="n", yaxt="n")
for(i in 1:100) {
	points(cos(2*i/100*pi), sin(2*i/100*pi),pch=20,col="gray")
}
 
mySample = sample(1:100,100)
for(i in 1:100) {
	found = FALSE
	nextItem = i
 
	# Pick a random color for this cycle
	color = sample(c(0:9,"A","B","C","D","E","F"),12,replace=T)
	lineColor = paste("#", paste(color[1:6],collapse=""),sep="")
 
	while(!found) {
		# Draw the cycle
 
		segments(cos(nextItem/50*pi), sin(nextItem/50*pi), cos(mySample[nextItem]/50*pi), sin(mySample[nextItem]/50*pi),col=lineColor,lwd=2)
		Sys.sleep(.4)
		if(mySample[nextItem] == i) {
			found = TRUE
		} else {
			nextItem = mySample[nextItem]
		}
	}
}

You can adjust the “Sys.sleep()” parameter to make the animation faster. I recommend running the code to see how the cycles “develop” over time, but here’s a snapshot of what I got: