Galton’s board all askew
Since Galton’s quincunx has fascinated me since the (early) days when I saw a model of it as a teenager in an industry museum near Birmingham, I jumped on the challenge to build an uneven nail version where the probabilities to end up in one of the boxes were not the Binomial ones. For instance, producing a uniform distribution with the maximum number of nails with probability ½ to turn right. And I obviously chose to try simulated annealing to figure out the probabilities, facing as usual the unpleasant task of setting the objective function, calibrating the moves and the temperature schedule. Plus, less usually, a choice of the space where the optimisation takes place, i.e., deciding on a common denominator for the (rational) probabilities. Should it be 2⁸?! Or more (since the solution with two levels also involves 1/3)? Using the functions
evol<-function(P){ Q=matrix(0,7,8) Q[1,1]=P[1,1];Q[1,2]=1-P[1,1] for (i in 2:7){ Q[i,1]=Q[i-1,1]*P[i,1] for (j in 2:i) Q[i,j]=Q[i-1,j-1]*(1-P[i,j-1])+Q[i-1,j]*P[i,j] Q[i,i+1]=Q[i-1,i]*(1-P[i,i]) Q[i,]=Q[i,]/sum(Q[i,])} return(Q)}
and
temper<-function(T=1e3){ bestar=tarP=targ(P<-matrix(1/2,7,7)) temp=.01 while (sum(abs(8*evol(R.01){ for (i in 2:7) R[i,sample(rep(1:i,2),1)]=sample(0:deno,1)/deno if (log(runif(1))/temp<tarP-(tarR<-targ(R))){P=R;tarP=tarR} for (i in 2:7) R[i,1:i]=(P[i,1:i]+P[i,i:1])/2 if (log(runif(1))/temp<tarP-(tarR<-targ(R))){P=R;tarP=tarR} if (runif(1)<1e-4) temp=temp+log(T)/T} return(P)}
I first tried running my simulated annealing code with a target function like
targ<-function(P)(1+.1*sum(!(2*P==1)))*sum(abs(8*evol(P)[7,]-1))
where P is the 7×7 lower triangular matrix of nail probabilities, all with a 2⁸ denominator, reaching
60
126 35
107 81 20
104 71 22 0
126 44 26 69 14
61 123 113 92 91 38
109 60 7 19 44 74 50
for 128P. With four entries close to 64, i.e. ½’s. Reducing the denominator to 16 produced once
8
12 1
13 11 3
16 7 6 2
14 13 16 15 0
15 15 2 7 7 4
8 0 8 9 8 16 8
as 16P, with five ½’s (8). But none of the solutions had exactly a uniform probability of 1/8 to reach all endpoints. Success (with exact 1/8’s and a denominator of 4) was met with the new target
(1+,1*sum(!(2*P==1)))*(.01+sum(!(8*evol(P)[7,]==1)))
imposing precisely 1/8 on the final line. With a solution with 11 ½’s
0.5
1.0 0.0
1.0 0.0 0.0
1.0 0.5 1.0 0.5
0.5 0.5 1.0 0.0 0.0
1.0 0.0 0.5 0.0 0.5 0.0
0.5 0.5 0.5 1.0 1.0 1.0 0.5
and another one with 12 ½’s:
0.5
1.0 0.0
1.0 .375 0.0
1.0 1.0 .625 0.5
0.5 0.5 0.5 0.5 0.0
1.0 0.0 0.5 0.5 0.0 0.5
0.5 1.0 0.5 0.0 1.0 0.5 0.0
Incidentally, Michael Proschan and my good friend Jeff Rosenthal have an 2009 American Statistician paper on another modification of the quincunx they call the uncunx! Playing a wee bit further with the annealing, and using a denominator of 840 let to a 60P with 13 ½’s out of 28
30
60 0
60 1 0
30 30 30 0
30 30 30 30 30
60 60 60 0 60 0
60 30 0 30 30 60 30
November 26, 2019 at 5:46 pm
Why no (Bayesian) two stage quincunx?
https://yihui.org/animation/example/quincunx2/
November 26, 2019 at 10:50 pm
Call it laziness, as I picked the first picture I found!!!
November 20, 2019 at 1:48 am
[…] November 20, 2019 By Donald Greer [This article was first published on R – Xi’an’s Og, and kindly contributed to R-bloggers]. (You can report issue about the content on this page […]
November 19, 2019 at 10:35 pm
[…] by data_admin [This article was first published on R – Xi’an’s Og, and kindly contributed to R-bloggers]. (You can report issue about the content on this page […]
November 19, 2019 at 3:56 pm
Yes! My introduction was the Eames exhibit in the IBM pavilion at the 1964 World’s Fair in NYC:
https://wikimili.com/en/Mathematica:_A_World_of_Numbers…_and_Beyond
I can watch them for hours.