## ABC-SAEM

Posted in Books, Statistics, University life with tags , , , , , , , , , , , , , , , , , , on October 8, 2019 by xi'an In connection with the recent PhD thesis defence of Juliette Chevallier, in which I took a somewhat virtual part for being physically in Warwick, I read a paper she wrote with Stéphanie Allassonnière on stochastic approximation versions of the EM algorithm. Computing the MAP estimator can be done via some adapted for simulated annealing versions of EM, possibly using MCMC as for instance in the Monolix software and its MCMC-SAEM algorithm. Where SA stands sometimes for stochastic approximation and sometimes for simulated annealing, originally developed by Gilles Celeux and Jean Diebolt, then reframed by Marc Lavielle and Eric Moulines [friends and coauthors]. With an MCMC step because the simulation of the latent variables involves an untractable normalising constant. (Contrary to this paper, Umberto Picchini and Adeline Samson proposed in 2015 a genuine ABC version of this approach, paper that I thought I missed—although I now remember discussing it with Adeline at JSM in Seattle—, ABC is used as a substitute for the conditional distribution of the latent variables given data and parameter. To be used as a substitute for the Q step of the (SA)EM algorithm. One more approximation step and one more simulation step and we would reach a form of ABC-Gibbs!) In this version, there are very few assumptions made on the approximation sequence, except that it converges with the iteration index to the true distribution (for a fixed observed sample) if convergence of ABC-SAEM is to happen. The paper takes as an illustrative sequence a collection of tempered versions of the true conditionals, but this is quite formal as I cannot fathom a feasible simulation from the tempered version and not from the untempered one. It is thus much more a version of tempered SAEM than truly connected with ABC (although a genuine ABC-EM version could be envisioned).

## Le Monde puzzle [#1092]

Posted in Statistics with tags , , , , , , , on April 18, 2019 by xi'an A Latin square Le Monde mathematical puzzle that I found rather dreary:

A hidden 3×3 board contains all numbers from 1 to 9. Anselm wants to guess the board and makes two proposals. Berenicke tells him how many entries are in the right rows and colums for each proposal, along with the information that no entry is at the right location. Anselm deduces the right board.

Which I solved by brute force and not even simulated annealing, first defining a target

ordoku1=ordoku2=matrix(1,9,2)
ordoku1[,1]=c(1,1,1,2,2,2,3,3,3)
ordoku1[,2]=rep(1:3,3)
ordoku2[,1]=c(3,2,3,1,2,3,2,1,1)
ordoku2[,2]=c(2,2,3,2,3,1,1,3,1)
fitz=function(ordo){
(sum(ordo[c(1,4,7),2]==1)==1)+(sum(ordo[c(2,5,8),2]==2)==1)+
(sum(ordo[c(3,6,9),2]==3)==0)+(sum(ordo[c(1,2,3),1]==1)==1)+
(sum(ordo[c(4,5,6),1]==2)==1)+(sum(ordo[c(7,8,9),1]==3)==2)+
(sum(ordo[c(6,7,9),2]==1)==2)+(sum(ordo[c(1,2,4),2]==2)==1)+
(sum(ordo[c(3,5,8),2]==3)==2)+(sum(ordo[c(4,8,9),1]==1)==1)+
(sum(ordo[c(7,2,5),1]==2)==1)+(sum(ordo[c(1,3,6),1]==3)==0)+
(!(0%in%apply((ordo-ordoku1)^2,1,sum)))+(!(0%in%apply((ordo-ordoku2)^2,1,sum)))
}


on a 9×9 board entry reproducing all items of information given by Berenicke. If all constraints are met, the function returns 14. And then searched for a solution at random:

temp=1
randw=function(ordo){
for (t in 1:1e6){
chlg=sample(1:9,2)
temp=ordo[chlg,]
ordo[chlg,]=ordo[chlg,]
ordo[chlg,]=temp
if (fitz(ordo)==14){
print(ordo);break()}}}


which produces the correct board

4 3 5
6 7 1
9 2 8


## Le Monde puzzle [#1088]

Posted in Books, Kids, R with tags , , , , , , , , on March 29, 2019 by xi'an A board (Ising!) Le Monde mathematical puzzle in the optimisation mode, again:

On a 7×7 board, what is the maximal number of locations that one can occupy when imposing at least two empty neighbours ?

Which I tried to solve by brute force and simulated annealing (what else?!), first defining a target

targ=function(tabz){
sum(tabz[-c(1,9),-c(1,9)]-1.2*(tabz[-c(1,9),-c(1,9)]*tabz[-c(8,9),-c(1,9)]
+tabz[-c(1,9),-c(1,9)]*tabz[-c(1,2),-c(1,9)]
+tabz[-c(1,9),-c(1,9)]*tabz[-c(1,9),-c(8,9)]
+tabz[-c(1,9),-c(1,9)]*tabz[-c(1,9),-c(1,2)]>2))}


on a 9×9 board where I penalise prohibited configuration by a factor 1.2 (a wee bit more than empty nodes). The perimeter of the 9×9 board is filled with ones and never actualised. (In the above convoluted products, the goal is to count how many neighbours of the entries equal to one are also equal to one. More than 2 is penalised.) The simulated annealing move is then updating the 9×9 grid gridz:

temp=1
maxarg=curarg=targ(gridz)
for (t in 1:1e3){
for (v in 1:1e4){
i=sample(2:8,1);j=sample(2:8,1)
newgrid=gridz;newgrid[i,j]=1-gridz[i,j]
newarg=targ(newgrid)
if (log(runif(1))<temp*(newarg-curarg)){
gridz=newgrid;curarg=newarg}}
temp=temp+.01}


and calls to the procedure always return 28 entries as the optimum, as in

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    1    0    1    0    1    0    1
[2,]    0    1    1    0    1    1    0
[3,]    1    1    0    1    0    1    1
[4,]    0    0    1    0    1    0    0
[5,]    1    1    0    1    0    1    1
[6,]    0    1    1    0    1    1    0
[7,]    1    0    1    0    1    0    1


As it happens, I had misread the wording of the original puzzle, which considered a dynamic placement of the units on the board, one at a time with two free neighbours imposed.

## Le Monde puzzle [#1087]

Posted in Books, Kids, R, Statistics with tags , , , , , on February 25, 2019 by xi'an

Given a (k,m) binary matrix, what is the maximum number S of entries with only one neighbour equal to one? Solve for k=m=2,…,13, and k=6,m=8.

For instance, for k=m=2, the matrix $\begin{matrix} 0 &0\\ 1 &1\\ \end{matrix}$

is producing the maximal number 4. I first attempted a brute force random filling of these matrices with only a few steps of explorations and got the numbers 4,8,16,34,44,57,… for the first cases. Since I was convinced that the square k² of a number k previously exhibited to reach its maximum as S=k² was again perfect in this regard, I then tried another approach based on Gibbs sampling and annealing (what else?):

gibzit=function(k,m,A=1e2,N=1e2){
temp=1 #temperature
board=sole=matrix(sample(c(0,1),(k+2)*(m+2),rep=TRUE),k+2,m+2)
board[1,]=board[k+2,]=board[,1]=board[,m+2]=0 #boundaries
maxol=counter(board,k,m) #how many one-neighbours?
for (t in 1:A){#annealing
for (r in 1:N){#basic gibbs steps
for (i in 2:(k+1))
for (j in 2:(m+1)){
prop=board
prop[i,j]=1-board[i,j]
u=runif(1)
if (log(u/(1-u))<temp*(counter(prop,k,m)-val)){
board=prop;val=counter(prop,k,m)
if (val>maxol){
maxol=val;sole=board}}
}}
temp=temp*2}
return(sole[-c(1,k+2),-c(1,m+2)])}


which leads systematically to the optimal solution, namely a perfect square k² when k is even and a perfect but one k²-1 when k is odd. When k=6, m=8, all entries can afford one neighbour exactly,

> gibzbbgiz(6,8)
 48
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    0    0    1    1    0    0    1
[2,]    1    0    0    0    0    0    0    1
[3,]    0    0    1    0    0    1    0    0
[4,]    0    0    1    0    0    1    0    0
[5,]    1    0    0    0    0    0    0    1
[6,]    1    0    0    1    1    0    0    1


but this does not seem feasible when k=6, m=7, which only achieves 40 entries with one single neighbour.

## Le Monde puzzle [#1075]

Posted in Books, Kids, R with tags , , , , , , , on November 22, 2018 by xi'an

A sequence of five integers can only be modified by subtracting an integer N from two neighbours of an entry and adding 2N to the entry.  Given the configuration below, what is the minimal number of steps to reach non-negative entries everywhere? Is this feasible for any configuration? As I quickly found a solution by hand in four steps, but missed the mathematical principle behind!, I was not very enthusiastic in trying a simulated annealing version by selecting the place to change inversely proportional to its value, but I eventually tried and also obtained the same solution:

      [,1] [,2] [,3] [,4] [,5]
-3    1    1    1    1
1   -1    1    1   -1
0    1    0    1   -1
-1    1    0    0    1
1    0    0    0    0


But (update!) Jean-Louis Fouley came up with one step less!

      [,1] [,2] [,3] [,4] [,5]
-3    1    1    1    1
3   -2    1    1   -2
2    0    0    1   -2
1    0    0    0    0


The second part of the question is more interesting, but again without a clear mathematical lead, I could only attempt a large number of configurations and check whether all admitted “solutions”. So far none failed.