Archive for simulated annealing

not a Bernoulli factory

Posted in Books, Kids, pictures, R with tags , , , , , , , on May 20, 2020 by xi'an

A Riddler riddle I possibly misunderstood:

Four isolated persons are given four fair coins, which can be either flipped once or returned without being flipped. If all flipped coins come up heads, the team wins! Else, if any comes up tails, or if no flip at all is done, it looses. Each person is further given an independent U(0,1) realisation. What is the best strategy?

Since the players are separated, I would presume the same procedure is used by all. Meaning that a coin is tossed with probability p, ie if the uniform is less than p, and untouched otherwise. The probability of winning is then

4(1-p)³p½+6(1-p)³p½²+4(1-p)p³½³+p⁴½⁴

which is maximum for p=0.3420391, with a winning probability of 0.2848424.

And an extra puzzle for free:

solve xxxx=2020

Where the integral part is the integer immediately below x. Puzzle that I first fail solving by brute force, because I did not look at negative x’s… Since the fourth root of 2020 is between 6 and 7, the solution is either x=6+ε or x=-7+ε, with ε in (0,1). The puzzle then becomes either

(6+ε)⌊(6+ε)⌊(6+ε)⌊6+ε⌋⌋ = (6+ε)⌊(6+ε)⌊36+6ε⌋⌋ = (6+ε)⌊(6+ε)(36+⌊6ε⌋)⌋ = 2020

where there are 6 possible integer values for ⌊6ε⌋, with only ⌊6ε⌋=5 being possible, turning the equation into

(6+ε)⌊41(6+ε) = (6+ε)(246+⌊41ε) = 2020

where again only ⌊42ε=40 being possible, ending up with

1716+286ε = 2020

which has no solution in (0,1). In the second case

(-7+ε)(-7+ε)(-7+ε)-7+ε⌋⌋ = (-7+ε)(-7+ε)(49+-7ε⌋)= 2020

shows that only -7ε=-3 is possible, leading to

(-7+ε)⌊46(-7+ε)) = (-7+ε) (-322+46ε⌋)=2020

with only 46ε=17 possible, hence

2135-305ε=2020

and

ε=115/305.

A brute force simulated annealing resolution returns x=-6.622706 after 10⁸ iterations. A more interesting question is to figure out the discontinuity points of the function

ℵ(x) = xxxx

as they seem to be numerous:

For instance, only 854 of the first 2020 integers enjoy a solution to ℵ(x)=n.

Le Monde puzzle [#1127]

Posted in Books, Kids, R, Statistics with tags , , , , on January 17, 2020 by xi'an

A permutation challenge as Le weekly Monde current mathematical puzzle:

When considering all games between 20 teams, of which 3 games have not yet been played, wins bring 3 points, losses 0 points, and draws 1 point (each). If the sum of all points over all teams and all games is 516, was is the largest possible number of teams with no draw in every game they played?

The run of a brute force R simulation of 187 purely random games did not produce enough acceptable tables in a reasonable time. So I instead considered that a sum of 516 over 187 games means solving 3a+2b=516 and a+b=187, leading to 142 3’s to allocate and 45 1’s. Meaning for instance this realisation of an acceptable table of game results

games=matrix(1,20,20);diag(games)=0
while(sum(games*t(games))!=374){
  games=matrix(1,20,20);diag(games)=0
  games[sample((1:20^2)[games==1],3)]=0}
games=games*t(games)
games[lower.tri(games)&games]=games[lower.tri(games)&games]*
    sample(c(rep(1,45),rep(3,142)))* #1's and 3'
    (1-2*(runif(20*19/2-3)<.5)) #sign
games[upper.tri(games)]=-games[lower.tri(games)]
games[games==-3]=0;games=abs(games)

Running 10⁶ random realisations of such matrices with no constraint whatsoever provided a solution with] 915,524 tables with no no-draws, 81,851 tables with 19 teams with some draws, 2592 tables with 18 with some draws and 3 tables with 17 with some draws. However, given that 10*9=90 it seems to me that the maximum number should be 10 by allocating all 90 draw points to the same 10 teams, and 143 3’s at random in the remaining games, and I reran a simulated annealing version (what else?!), reaching a maximum 6 teams with no draws. Nowhere near 10, though!

riddle on a circle

Posted in Books, Kids, R, Travel with tags , , , , , , , on December 22, 2019 by xi'an

The Riddler’s riddle this week provides another opportunity to resort to brute-force simulated annealing!

Given a Markov chain defined on the torus {1,2,…,100} with only moves a drift to the right (modulo 100) and a uniformely random jump, find the optimal transition matrix to reach 42 in a minimum (average) number of moves.

Which I coded in my plane to Seattle, under the assumption that there is nothing to do when the chain is already in 42. And the reasoning that there is not gain (on average) in keeping the choice between right shift and random jump random.

dure=min(c(41:0,99:42),50)
temp=.01
for (t in 1:1e6){
  i=sample((1:100)[-42],1)
  dura=1+mean(dure)
  if (temp*log(runif(1))<dure[i]-dura) dure[i]=dura
  if(temp*log(runif(1))<dure[i]-(dura<-1+dure[i*(i<100)+1])) 
    dure[i]=dura 
  temp=temp/(1+.1e-4*(runif(1)>.99))}

In all instances, the solution is to move at random for any position but those between 29 and 41, for an average 13.64286 number of steps to reach 42. (For values outside the range 29-42.)

Galton’s board all askew

Posted in Books, Kids, R with tags , , , , , , , on November 19, 2019 by xi'an

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

stochastic magnetic bits, simulated annealing and Gibbs sampling

Posted in Statistics with tags , , , , , , , , , on October 17, 2019 by xi'an

A paper by Borders et al. in the 19 September issue of Nature offers an interesting mix of computing and electronics and optimisation. With two preparatory tribunes! One [rather overdone] on Feynman’s quest. As a possible alternative to quantum computers for creating probabilistic bits. And making machine learning (as an optimisation program) more efficient. And another one explaining more clearly what is in the paper. As well as the practical advantage of the approach over quantum computing. As for the paper itself, the part I understood about factorising an integer F via minimising the squared difference between a product of two integers and F and using simulated annealing sounded rather easy, while the part I did not about constructing a semi-conductor implementing this stochastic search sounded too technical (especially in the métro during rush hour). Even after checking the on-line supplementary material. Interestingly, the paper claims for higher efficiency thanks to asynchronicity than a regular Gibbs simulation of Boltzman machines, quoting Roberts and Sahu (1997) without further explanation and possibly out of context (as the latter is not concerned with optimisation).