**T**oday I received in the mail a copy of the short book published by edp sciences after the courses we gave last year at the astrophysics summer school, in Autrans. Which contains a quick introduction to ABC extracted from my notes (which I still hope to turn into a book!). As well as a longer coverage of Bayesian foundations and computations by David Stenning and David van Dyk.

## Archive for the R Category

## ABC intro for Astrophysics

Posted in Books, Kids, Mountains, R, Running, Statistics, University life with tags ABC, Approximate Bayesian computation, Autrans, Bayesian foundations, Bayesian methodology, Book, computational astrophysics, review, Statistics for Astrophysics, summer course, survey, Vercors on October 15, 2018 by xi'an## Le Monde puzzle [#1070]

Posted in Books, Kids, R, University life with tags CEREMADE, competition, Dauphine, dynamic programming, Le Monde, mathematical puzzle, optimisation, R on October 11, 2018 by xi'an**R**ewording Le Monde mathematical puzzle fifth competition problem

For the 3×3 tables below, what are the minimal number of steps to move from left to rights when the yellow tokens can only be move to an empty location surrounded by two other tokens?

In the 4×4 table below, there are 6 green tokens. How many steps from left to right?

Painful and moderately mathematical, once more… For the first question, a brute force simulation of random valid moves of length less than 100 returns solutions in 4 steps:

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

But this is not an acceptable move because of the “other” constraint. Imposing this constraint leads to a solution in 9 steps, but is this the lowest bound?! It actually took me most of the weekend (apart from a long drive to and from a short half-marathon!) to figure out a better strategy than brute force random exploration: the trick I eventually figured out is to start from the finishing (rightmost) value F of the grid and look at values with solutions available in 1,2,… steps. This is not exactly dynamic programming, but it keeps the running time under control if there is a solution associated with the starting (leftmost) value S. (Robin proceeded reverse-wise, which in retrospect is presumably equivalent, if faster!) The 3×3 grid has 9 choose 5, ie 126, possible configurations with 5 coins, which means the number of cases remains under control. And even so for the 4×4 grid with 6 coins, with 8008 configurations. This led to a 9 step solution for n=3 and the proposed starting grid in yellow:

[1] 1 1 1 0 0 1 0 0 1 [1] 1 1 0 0 1 1 0 0 1 [1] 1 1 0 1 1 0 0 0 1 [1] 0 1 0 1 1 1 0 0 1 [1] 0 1 1 1 0 1 0 0 1 [1] 1 1 1 1 0 0 0 0 1 [1] 0 1 1 1 1 0 0 0 1 [1] 0 0 1 1 1 0 0 1 1 [1] 0 0 1 1 0 0 1 1 1 [1] 0 0 1 0 0 1 1 1 1

and a 19 step solution for n=4:

[1] 1 0 0 0 0 1 0 0 0 1 1 1 0 0 0 1 [1] 1 0 0 0 0 1 0 0 0 0 1 1 0 0 1 1 [1] 1 0 0 0 0 1 1 0 0 0 1 1 0 0 1 0 [1] 1 1 0 0 0 1 1 0 0 0 0 1 0 0 1 0 [1] 1 1 0 0 0 0 1 1 0 0 0 1 0 0 1 0 [1] 1 1 1 0 0 0 1 1 0 0 0 0 0 0 1 0 [1] 1 0 1 1 0 0 1 1 0 0 0 0 0 0 1 0 [1] 1 1 1 1 0 0 1 0 0 0 0 0 0 0 1 0 [1] 1 1 0 1 0 1 1 0 0 0 0 0 0 0 1 0 [1] 1 0 0 1 1 1 1 0 0 0 0 0 0 0 1 0 [1] 0 0 0 1 1 1 1 0 0 0 1 0 0 0 1 0 [1] 0 0 0 1 1 1 0 0 0 1 1 0 0 0 1 0 [1] 0 0 0 1 1 1 0 0 1 1 0 0 0 0 1 0 [1] 0 0 0 1 0 1 0 0 1 1 0 0 0 1 1 0 [1] 0 0 0 1 0 1 0 0 1 0 0 0 1 1 1 0 [1] 0 0 0 1 1 1 0 0 1 0 0 0 1 0 1 0 [1] 0 0 0 1 0 1 0 0 1 1 0 0 1 0 1 0 [1] 0 0 0 1 0 1 0 0 0 1 1 0 1 0 1 0 [1] 0 0 0 1 0 1 1 0 0 1 1 0 1 0 0 0

The first resolution takes less than a minute and the second one just a few minutes (or less than my short morning run!). Surprisingly, using this approach does not require more work, which makes me wonder at the solution Le Monde journalists will propose. Given the (misguided) effort put into my resolution, seeing a larger number of points for this puzzle is no wonder.

## Le Monde puzzle [#1068]

Posted in Books, Kids, R with tags competition, grid, Le Monde, mathematical puzzle, R, simulated annealing, stochastic algorithms on October 3, 2018 by xi'an**A **purely (?) algorithmic Le Monde mathematical puzzle

For the table below, what is the minimal number of steps required to reach equal entries when each step consists in adding ones to three entries sitting in aL, such as (7,11,12) or (5,6,10)? Same question for the inner table of four in yellow.

For the inner table, this is straightforward as there are four possible L’s, three equations like 6+n⁶=7+n⁷, and two degrees of freedom leading to a unique entry of N=13 (impossible!) or 16 (feasible). Hence adding 10 L’s. For the entire table, summing up all entries after completion leads to 16N, which is also equal to 1+3+3+…+16+M, where M is the number of added L’s, itself equal to 138+3O, if O denotes the number of ones added. Hence M can only take the values 18, 21, … It took me quite a while to R code an approach to complete the table into 16 18’s, as my versions of simulated annealing did not seem to converge. In the end, I used a simplified version where the table was completed by multinomial draws, like a M(17;3⁻¹,3⁻¹,3⁻¹) for the upper left corner, corresponding to random draws of one of the 36 available L’s, which should be used 50 times in total, and then augmented or reduced of one L depending on the value at a randomly selected entry. Leading to the result

> aneal(grid=c(1,3,3:13,15,15,16),maxT=1e5) [1] 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18

## Le Monde puzzle [#1068]

Posted in Books, Kids, R with tags CEREMADE, competition, Dauphine, Le Monde, mathematical puzzle, optimisation, R, simulated annealing on September 26, 2018 by xi'an**A**nd here is the third Le Monde mathematical puzzle open for competition:

Consider for this puzzle only integers with no zero digits. Among these an integer x=a¹a²a³… isrefinedif it is a multiple of its scion, the integer defined as x without the first digit, y=a²a³…. Find the largest refined integer x such the sequence of the successive scions of x with more than one digit is entirely refined. An integer x=a¹a²a… is distilled if it is a multiple of its grand-scion, the integer defined as x without the first two digits, z=a³… Find the largest distilled integer x such the sequence of the successive scions of x with more than two digits is entirely distilled.

Another puzzle amenable to a R resolution by low-tech exploration of possible integers, first by finding refined integers, with no solution between 10⁶ and 10⁸ [*meaning there is no refined integer larger than 10⁶*] and then checking which refined integers have refined descendants, e.g.,

raf=NULL for (x in (1e1+1):(1e6-1)){ y=x%%(10^trunc(log(x,10))) if (y>0){ if (x%%y==0) raf=c(raf,x)}}

that returns 95 refined integers. And then

for (i in length(raf):1){ gason=raf[i] keep=(gason%in%raf) while (keep&(gason>100)){ gason=gason%%(10^trunc(log(gason,10))) keep=keep&(gason%in%raf)} if (keep) break()}}

that returns 95,625 as the largest refined integer with the right descendance. Rather than finding all refined integers at once, going one digit at a time and checking that some solutions have proper descendants is actually faster.

Similarly, running an exploration code up to 10⁹ produces 1748 distilled integers, with maximum 9,977,34,375, so it is unlikely this is the right upper bound but among these the maximum with the right distilled descendants is 81,421,875. As derived by

rad=(1:99)[(1:99)%%10>0] for (dig in 2:12){ for (x in ((10^dig+1):(10^{dig+1}-1))){ y=x%%(10^{dig-1}) if (y>0){ if (x%%y==0){ if (min(as.integer(substring(x,seq(nchar(x)),seq(nchar(x)))))>0){ rad=c(rad,x) y=x%%(10^dig) keep=(y%in%rad) while (keep&(y>1e3)){ y=y%%(10^trunc(log(y,10))) keep=keep&(y%in%rad)} if (keep) solz=x}}}} if (solz<10^dig) break() topsol=max(solz)}

## Le Monde puzzle [#1067]

Posted in Books, Kids, R with tags Al-Kashi's Theorem, Alice and Bob, competition, geometry, law of cosines, Le Monde, mathematical puzzle, trigonometry on September 19, 2018 by xi'an**T**he second Le Monde mathematical puzzle in the new competition is sheer trigonometry:

When in the above figures both triangles ABC are isosceles and the brown segments are all of length 25cm, find the angle in A and the value of DC², respectively.

This could have been solved by R coding the various possible angles of the three segments beyond BC until the isosceles property is met, but it went much faster by pen and paper, leading to an angle of π/9 in the first case and a square of 1250 in the second case. The third puzzle is basic arithmetic that only seems solvable by enumeration…

## Le Monde puzzle [#1066]

Posted in Books, Kids, R with tags CEREMADE, competition, Dauphine, Le Monde, mathematical puzzle, optimisation, R, simulated annealing on September 13, 2018 by xi'an**R**ecalling Le Monde mathematical puzzle first competition problem

For the X table below, what are the minimal number of lights that are on (green) to reach the minimal and maximal possible numbers of entries (P) with an even (P as pair) number of neighbours with lights on? In the illustration below, there are 16 lights on (green) and 31 entries with an even number of on-neighbours.

*As suggested last week, this was amenable to a R resolution by low-tech simulated annealing although the number of configurations was not that large when accounting for symmetries. The R code is a wee bit long for full reproduction here but it works on moving away from a random filling of this cross by 0’s and 1’s, toward minimising or maximising the number of P’s, this simulated annealing loop being inserted in another loop recording the minimal number of 1’s in both cases. A first round produced 1 and 44 for the minimal and maximal numbers of P’s, respectively, associated with at least 16 and 3 1’s, respectively, but a longer run exhibited 45 for 6 1’s crossing one of the diagonals of the X, made of two aligned diagonals of the outer 3×3 tables. (This was confirmed by both Amic and Robin in Dauphine!) The next [trigonometry] puzzle is on!*

## riddles on a line [#2]

Posted in Books, Kids, R with tags British Columbia, FiveThirtyEight, mathematical puzzle, probability, R, The Riddler, Tofino, Vancouver Island on September 11, 2018 by xi'an**A** second Riddle(r), with a puzzle related with the integer set Ð={,12,3,…,N}, in that it summarises as

Given a random walk on Ð, starting at the middle N/2, with both end states being absorbing states, and a uniform random move left or right of the current value to the (integer) middle of the corresponding (left or right) integer interval, what is the average time to one absorbing state as a function of N?

Once the Markov transition matrix M associated with this random walk is defined, the probability of reaching an absorbing state in t steps can be derived from the successive powers of M by looking at the difference between the probabilities to be (already) absorbed at both t-1 and t steps. From which the average can be derived.

avexit <- function(N=100){ #transition matrix M for the walk #1 and N+2 are trapping states tranz=matrix(0,N+2,N+2) tranz[1,1]=tranz[N+2,N+2]=1 for (i in 2:(N+1)) tranz[i,i+max(trunc((N+1-i)/2),1)]=tranz[i,i-max(trunc((i-2)/2),1)]=1/2 #probabilities of absorption prowin=proloz=as.vector(0) init=rep(0,N+2) init[trunc((N+1)/2)]=1 #first position curt=init while(1-prowin[length(prowin)]-proloz[length(prowin)]>1e-10){ curt=curt%*%tranz prowin=c(prowin,curt[1]) proloz=c(proloz,curt[N+2])} #probability of new arrival in trapping state probz=diff(prowin+proloz) return(sum((2:length(proloz))*probz))}

leading to an almost linear connection between N and expected trapping time.