## splitting a field by annealing

Posted in Kids, pictures, R, Statistics with tags , , , , , , , , on October 18, 2017 by xi'an

A recent riddle [from The Riddle] that I pondered about during a [long!] drive to Luxembourg last weekend was about splitting a square field into three lots of identical surface for a minimal length of separating wire… While this led me to conclude that the best solution was a T like separation, I ran a simulated annealing R code on my train trip to AutransValence, seemingly in agreement with this conclusion.I discretised the square into n² units and explored configurations by switching two units with different colours, according to a simulated annealing pattern (although unable to impose connectivity on the three regions!):

partz=matrix(1,n,n)
partz[,1:(n/3)]=2;partz[((n/2)+1):n,((n/3)+1):n]=3
#counting adjacent units of same colour
nood=hood=matrix(4,n,n)
for (v in 1:n2) hood[v]=bourz(v,partz)
minz=el=sum(4-hood)
for (t in 1:T){
colz=sample(1:3,2) #picks colours
a=sample((1:n2)[(partz==colz[1])&(hood<4)],1)
b=sample((1:n2)[(partz==colz[2])&(hood<4)],1)
partt=partz;partt[b]=colz[1];partt[a]=colz[2]
#collection of squares impacted by switch
nood=hood
voiz=unique(c(a,a-1,a+1,a+n,a-n,b-1,b,b+1,b+n,b-n))
voiz=voiz[(voiz>0)&(voiz<n2)]
for (v in voiz) nood[v]=bourz(v,partt)
if (nood[a]*nood[b]>0){
difz=sum(nood)-sum(hood)
if (log(runif(1))<difz^3/(n^3)*(1+log(10*rep*t)^3)){
el=el-difz;partz=partt;hood=nood
if (el<minz){ minz=el;cool=partz}
}}}


(where bourz computes the number of neighbours), which produces completely random patterns at high temperatures (low t) and which returns to the T configuration (more or less):if not always, as shown below:Once the (a?) solution was posted on The Riddler, it appeared that one triangular (Y) version proved better than the T one [if not started from corners], with a gain of 3% and that a curved separation was even better with an extra gain less than 1% [solution that I find quite surprising as straight lines should improve upon curved ones…]

## sequence riddle

Posted in Kids, R with tags , , , , , on August 10, 2017 by xi'an

The riddle this week on The Riddler was about finding the largest sequence of integers between 1 and 100 such that each integer is only used once and always followed by a multiple or a factor. A basic R code searching at random [and programmed during a massive downpour on Skye] led to a solution of 69:

although there is no certainty this is the best p… And the solutions posted the next week showed sequences with length 77! [Interestingly, both posted solutions have a sequence starting with 87. And they seem to exploit the graph of connections between integers in a much more subtle way that my random exploration of subsequences.]

## easy riddle

Posted in Books, Kids, R with tags , , , , , on July 12, 2017 by xi'an

From the current Riddler, a problem that only requires a few lines of code and a few seconds of reasoning. Or not.

N households each stole the earnings from one of the (N-1) other households, one at a time. What is the probability that a given household is not burglarised? And what are the expected final earnings of each household in the list, assuming they all start with \$1?

The first question is close to Feller’s enveloppe problem in that

$\left(1-\frac{1}{N-1}\right)^{N-1}$

is close to exp(-1) for N large. The second question can easily be solved by an R code like

N=1e3;M=1e6
fina=rep(1,N)
for (v in 1:M){
ordre=sample(1:N)
vole=sample(1:N,N,rep=TRUE)
while (min(abs(vole-(1:N)))==0)
vole[abs(vole-(1:N))==0]=sample(1:N,
sum(vole-(1:N)==0))
cash=rep(1,N)
for (t in 1:N){
cash[ordre[t]]=cash[ordre[t]]+cash[vole[t]];cash[vole[t]]=0}
fina=fina+cash[ordre]}


which returns a pretty regular exponential-like curve, although I cannot figure the exact curve beyond the third burglary. The published solution gives the curve

${\frac{N-2}{N-1}}^{999}\times 2+{\frac{1}{N-1}}^{t-1}\times{\frac{N-1}{N}}^{N-t}\times\frac{N}{N-1}$

corresponding to the probability of never being robbed (and getting on average an extra unit from the robbery) and of being robbed only before robbing someone else (with average wealth N/(N-1)).

## [un]solved riddles

Posted in Books, Kids, R with tags , , on July 4, 2017 by xi'an

On the Riddler of last week, first a birthday puzzle:

Given a group of 23 persons, what is the probability of observing three pairs of identical birthdays?

which can be found by a quick simulation as

ave=0
for (t in 1:1e6){
dupz=dates[duplicated(sample(1:365,23,rep=TRUE))]
ave=ave+as.integer((length(dupz)==3)&
(length(unique(dupz))==3))}}
ave/M


returning a value of 0.0183, but which combinatoric resolution I could not fully fathom without a little help from a friend (-ly blog). I had the multinomial coefficient

${23\choose 2\ 2\ 2}$

for the allocation of the 23 persons to one of the three pairs or none, as well as the probability

$\dfrac{365\times\cdots\times(365-(23-4)+1)}{365^{23}}$

but I had forgotten the 3! in the denominator for the permutations of the three pairs, which leads again to

${23 \choose 2\ 2\ 2}\dfrac{365\times\cdots\times(365-(23-4)+1)}{3!\times365^{23}} = 0.0183$

A question that also led to an unsolved question: in the even this probability was much smaller, what is an easy way into constructing a more efficient importance sampler?

The second riddle was just as easy to code in R:

A game of tag goes by the following rules: (i) anyone untagged can tag anyone untagged; (ii) anyone tagged by a player tagged gets untagged; (iii) the winner is the last untagged player. What is the expected number of runs for N players?

The outcome of

game=function(N=12){
a=rep(0,N);T=0
while (sum(a==0)>1){
ij=sample((1:N)[a==0],2)
a[ij[2]]=ij[1];a[a==ij[2]]=0
T=T+1}
return(T)}


leads to an average value of

$2^{N-1}-1$

but I had no clear quick explanation for the doubling phenomenon. Until I picked a pen and a sheet of paper and drew the last steps of the game: to decrease down to 1, the size of the untagged survivors has to get through …,3,2 and each time the eliminated player needs to have tagged no other player since otherwise the population grows again. This has to apply all the way to the second round, where N-1 players remain and the one tagged needs to be anyone but the one who tagged the first one. And so on…

## non-Bayesian decision riddle

Posted in Statistics with tags , , , on June 22, 2017 by xi'an

As a continuation of the Bayesian resolution of last week riddle, I looked at a numeric resolution of the four secretaries problem, while in the train back from Rouen (and trying to block the chatter of my neighbours, a nuisance I find myself more and more sensitive to!). The target function is defined as

gainz=function(b,c,T=1e4,type="raw"){
x=matrix(runif(4*T),ncol=4)
maz=t(apply(x,1,cummax))
zam=t(apply(x[,4:1],1,cummax))
if (type=="raw"){return(mean(
((x[,2]>b*x[,1])*x[,2]+
(x[,2]<b*x[,1])*((x[,3]>c*maz[,2])*x[,3]+
(x[,3]<c*maz[,2])*x[,4]))/maz[,4]))}
if (type=="global"){return(mean(
((x[,2]>b*x[,1])*(x[,2]==maz[,4])+
(x[,2]<b*x[,1])*((x[,3]>c*maz[,2])*(x[,3]==maz[,4])+
(x[,3]<c*maz[,2])*(x[,4]==maz[,4])))))}
if (type=="remain"){return(mean(
((x[,2]>b*x[,1])*(x[,2]==zam[,3])+
(x[,2]<b*x[,1])*((x[,3]>c*maz[,2])*(x[,3]==zam[,2])+
(x[,3]<c*maz[,2])*(x[,4]==zam[,2])))))}}


where the data is generated from a U(0,1) distribution as the loss functions are made scaled free by deciding to always sacrifice the first draw, x¹. This function is to be optimised in (b,c) and hence I used a plain vanilla simulated annealing R code:

avemale=function(T=3e4,type){
b=c=.5
maxtar=targe=gainz(b,c,T=1e4,type)
temp=0.1
for (t in 1:T){
bp=b+runif(1,-temp,temp)
cp=c+runif(1,-temp,temp)
parge=(bp>0)*(cp>0)*gainz(bp,cp,T=1e4,type)
if (parge>maxtar){
b=bs=bp;c=cs=cp;maxtar=targe=parge}else{
if (runif(1)<exp((parge-targe)/temp)){
b=bp;c=cp;targe=parge}}
temp=.9999*temp}
return(list(bs=bs,cs=cs,max=maxtar))}


with outcomes

• b=1, c=.5, and optimum 0.8 for the raw type
• b=c=1 and optimum 0.45 for the global type
• b undefined, c=2/3 and optimum 0.75 for the remain type

## Bayesian decision riddle

Posted in Books, Kids, Statistics with tags , , , , on June 15, 2017 by xi'an

The current puzzle on The Riddler is a version of the secretary problem with an interesting (?) Bayesian solution.

Given four positive numbers x¹, x², x³, x⁴, observed sequentially, the associated utility is the value of x at the stopping time. What is the optimal stopping rule?

While nothing is mentioned about the distribution of the x’s, I made the assumption that they were iid and uniformly distributed over (0,M), with M unknown and tried a Bayesian resolution with the non-informative prior π(M)=1/M. And failed. The reason for this failure is that the expected utility is infinite at the first step: while the posterior expected utility is finite with three and two observations, meaning I can compare stopping and continuing at the second and third steps, the predicted expected reward for continuing after observing x¹ does not exist because the expected value of max(x¹,x²) given x¹ does not exist. As the predictive density of x² is max(x¹,x²)⁻²…  Several alternatives are possible to bypass this impossible resolution, from changing the utility function to picking another reference prior.

For instance, using a prior like π(M)=1/M² l(and the same monetary return utility) leads to a proper optimal solution, namely

1. always wait for the second observation x²
2. stop at x² if x²>11x¹/12, else wait for x³
3. stop at x³ if x³>23 max(x¹,x²)/24, else observe x⁴

obtained analytically on a bar table in Rouen (and checked numerically later).

Another approach is to try to optimise the probability to pick the largest amount of the four x’s, but this is not leading to an interesting solution, since it corresponds to picking the first maximum after x¹, while picking the largest among remaining ones leads to a somewhat convoluted solution I have no patience to produce here! Plus this is not a really pertinent loss function as it does not discriminate enough against waiting…

## continental divide

Posted in Books, Kids, pictures, R with tags , , , , on May 19, 2017 by xi'an

While the Riddler puzzle this week was anticlimactic,  as it meant filling all digits in the above division towards a null remainder, it came as an interesting illustration of how different division is taught in the US versus France: when I saw the picture above, I had to go and check an American primary school on-line introduction to division, since the way I was taught in France is something like that

with the solution being that 12128316 = 124 x 97809… Solved by a dumb R exploration of all constraints:

for (y in 111:143)
for (z4 in 8:9)
for (oz in 0:999){
z=oz+7e3+z4*1e4
x=y*z
digx=digits(x)
digz=digits(z)
if ((digz[2]==0)&(x>=1e7)&(x<1e8)){
r1=trunc(x/1e4)-digz[5]*y
if ((digz[5]*y>=1e3)&(digz[4]*y<1e4) &(r1>9)&(r1<100)){
r2=10*r1+digx[4]-7*y
if ((7*y>=1e2)&(7*y<1e3)&(r2>=1e2)&(r2<1e3)){
r3=10*r2+digx[3]-digz[3]*y
if ((digz[3]*y>=1e2)&(digz[3]*y<1e3)&(r3>9)&(r3<1e2)){
r4=10*r3+digx[2]
if (r4<y) solz=rbind(solz,c(y,z,x))
}}}}


Looking for a computer-free resolution, the constraints on z exhibited by the picture are that (a) the second digit is 0 and the fourth digit is 7.  Moreover, the first and fifth digits are larger than 7 since y times these digits is a four-digit number. Better, since the second subtraction from a three-digit number by 7y returns a three-digit number and the third subtraction from a four-digit number by ny returns a two-digit number, n is larger than 7 but less than the first and fifth digits. Ergo, z is necessarily 97809! Furthermore, 8y<10³ and 9y≥10³, which means 111<y<125. Plus the constraint that 1000-8y≤99 implies y≥112. Nothing gained there! This leaves 12 values of y to study, unless there is another restriction I missed…