Jun 10

Repulsive dots pattern, the difference of distance

What if you wanted to randomly place objects into a field, and the more objects you had, the more they rejected newcomers placed nearby? To find out, I setup a simulation. The code, shown at the end, isn’t all that interesting, and the plots shown below aren’t all that special. I think there is one interesting part of this, and that’s how the clustering changes depending on how distance is measured. One of the plots uses the traditional “L2″ distance, the other uses L1” (Manhattan taxi cab) measure . Each plot shown below has almost exactly the same number of dots (277 vs 279). Can you tell which uses L1 and which uses L2 just by looking?

Plot A:

Plot B:

Here’s the code. Run it and see for yourself. Make sure to change adjust the values which have comments next to them. Uncommenting “print(force)” can help you pack a maxRepulse value.

calcRepulse <- function(x,y,dots,use="L2") {
	force = 0
	i = 1
	while(i <= dim(dots)[1] && dots[i,1] != 0) {
		if(use == "L2") {
			force = force + 1/( (x-dots[i,1])^2 + (y-dots[i,2])^2 )
		} else if(use == "L1") {
			force = force + 1/( abs(x-dots[i,1]) + abs(y-dots[i,2]) )
		i = i+1
	# print(force)
place = 1 #Maximum number of dots to place, change this to something bigger
dots = matrix(rep(0,place*2),ncol=2)
maxTries = place * 10
maxRepulse = 1 # Anything above this will be rejected as too repulsive a location
dist2use = "" # Pick L1 or L2
placed = 0
tries = 0
while(placed < place && tries < maxTries) {
	x = runif(1)
	y = runif(1)
	if(calcRepulse(x,y,dots,dist2use) < maxRepulse) {
		dots[(placed + 1),1] = x
		dots[(placed + 1),2] = y
		placed = placed + 1
	tries = tries + 1

May 10

R: More plotting fun with Poission

Coded as follows:

x = seq(.001,50,.001)

May 10

First annual R plot replication prize

Click image for the full-size version

$100 to the first person who can figure out how I created this plot and replicate it. Some hints:

  • It was done in R.
  • There is only one underlying probability distribution involved (one “rdist()“).
  • Including the “plot” statement, I created this with 3 short lines of code.

This is based on a random sampling of unstated size, so I don’t expect that your graph will be an absolute, exact match.

I’ll add $1 to the prize for every day that goes by without a winner until the end of the year. After that I’ll consider it an unsolved mystery and reveal the code I used.

Post your guesses for the code as a comment to this post. First correct answer wins. Good luck to all!

Apr 10

R: more plotting fun, this time with the Poisson

Click on image for a larger version. Here is the code:

par(mar=c(0,0,0,0)) plot(sort(rpois(10000,100))/rpois(10000,100),frame.plot=F,pch=20,col="blue")

Apr 10

R: another nifty graph

Make sure to click on the image to see the large version. Code for this graph:

moxbuller = function(n) {   
	u = runif(n)   
	v = runif(n)   
	x = cos(2*pi*u)*sqrt(-2*log(v))  
	y = sin(2*pi*v)*sqrt(-2*log(u))
	r = list(x=x, y=y)
r = moxbuller(50000) 
plot(r$x,r$y, pch=".", col="blue", cex=1.2)