**A**s the Le Monde mathematical puzzle of this week was a geometric one *(the quadrangle ABCD is divided into two parts with the same area, &tc…)* , with no clear R resolution, I chose to bypass it. In this April 3 issue, several items of interest: first, a report by Etienne Ghys on Yakov Sinaï’s Abel Prize for his work “between determinism and randomness”, centred on ergodic theory for dynamic systems, which sounded like the ultimate paradox the first time I heard my former colleague Denis Bosq give a talk about it in Paris 6. Then a frightening fact: the summer conditions have been so unusually harsh in Antarctica (or at least near the Dumont d’Urville French austral station) that none of the 15,000 Adélie penguin couples studied there managed to keep their chick alive. This was due to an ice shelf that did not melt at all over the summer, forcing the penguins to walk an extra 40k to reach the sea… Another entry on the legal obligation for all French universities to offer a second chance exam, no matter how students are evaluated in the first round. (Too bad, I always find writing a second round exam a nuisance.)

## Archive for R

## Le Monde sans puzzle [& sans penguins]

Posted in Books, Kids, R, University life with tags intToBits(), Le Monde, mathematical puzzle, R, StackExchange, stackoverflow on April 12, 2014 by xi'an## Le Monde puzzle [#860]

Posted in Books, Kids, R with tags awalé, Le Monde, mathematical puzzle, R, recursive function on April 4, 2014 by xi'an**A** Le Monde mathematical puzzle that connects to my awalé post of last year:

For N≤18, N balls are placed in N consecutive holes. Two players, Alice and Bob, consecutively take two balls at a time provided those balls are in contiguous holes. The loser is left with orphaned balls. What is the values of N such that Bob can win, no matter what is Alice’s strategy?

**I** solved this puzzle by the following R code that works recursively on N by eliminating all possible adjacent pairs of balls and checking whether or not there is a winning strategy for the other player.

topA=function(awale){ # return 1 if current player can win, 0 otherwise best=0 if (max(awale[-1]*awale[-N])==1){ #there are adjacent balls remaining for (i in (1:(N-1))[awale[1:(N-1)]==1]){ if (awale[i+1]==1){ bwale=awale bwale[c(i,i+1)]=0 best=max(best,1-topA(bwale)) } }} return(best) } for (N in 2:18) print(topA(rep(1,N)))

which returns the solution

[1] 1 [1] 1 [1] 1 [1] 0 [1] 1 [1] 1 [1] 1 [1] 0 [1] 1 [1] 1 [1] 1 [1] 1 [1] 1 [1] 0 [1] 1 [1] 1 [1] 1 <pre>

(brute-force) answering the question that N=5,9,15 are the values where Alice has no winning strategy if Bob plays in an optimal manner**.** (The case N=5 is obvious as there always remains two adjacent 1′s once Alice removed any adjacent pair. The case N=9 can also be shown to be a lost cause by enumeration of Alice’s options.)

## Le Monde puzzle [#857]

Posted in Books, Kids, R with tags is.prim, Le Monde, mathematical puzzle, number theory, prime factor decomposition, prime.factor, R, schoolmath on March 22, 2014 by xi'an**A** rather bland case of Le Monde mathematical puzzle :

Two positive integers x and y are turned into s=x+y and p=xy. If Sarah and Primrose are given S and P, respectively, how can the following dialogue happen?

I am sure you cannot find my numberNow you told me that, I can, it is 46.

and what are the values of x and y?

**I**n the original version, it was unclear whether or not each person knew she had the sum or the product. Anyway, the first person in the dialogue has to be Sarah, since a product p equal to a prime integer would lead Primrose to figure out x=1 and hence s=p+1. (Conversely, having observed the sum s cannot lead to deduce x and y.) This means x+y-1 is *not* a prime integer. Now the deduction of Primrose that the sum is 46 implies p can be decomposed only once in a product such that x+y-1 is not a prime integer. If p=45, this is the case since 45=15×3 and 45=5×9 lead to 15+3-1=17 and 5+9-1=13, while 45=45×1 leads to 45+1-1=45. Other solutions fail, as demonstrated by the R code:

> for (x in 1:23){ + fact=c(1,prime.factor(x*(46-x))) + u=0; + for (i in 1:(length(fact)-1)) + u=u+1-is.prim(prod(fact[1:i])+prod(fact[-(1:i)])-1) + if (u==1) print(x)} [1] 1

**B**usser and Cohen argue much more wisely in their solution that any non-prime product p other than 45 would lead to p+1 as an acceptable sum s, hence would prevent Primrose from guessing s.

## fine-sliced Poisson [a.k.a. sashimi]

Posted in Books, Kids, pictures, R, Running, Statistics, University life with tags convergence assessment, ENSAE, Gibbs sampling, Monte Carlo Statistical Methods, Poisson distribution, R, sample, sampling from an atomic population, slice sampling, slow convergence on March 20, 2014 by xi'an**A**s my student Kévin Guimard had not mailed me his own Poisson slice sampler of a Poisson distribution, I could not tell why the code was not working! My earlier post prompted him to do so and a somewhat optimised version is given below:

nsim = 10^4 lambda = 6 max.factorial = function(x,u){ k = x parf=1 while (parf*u<1){ k = k + 1 parf = parf * k } k = k - (parf*u>1) return (k) } x = rep(floor(lambda), nsim) for (t in 2:nsim){ v1 = ceiling((log(runif(1))/log(lambda))+x[t-1]) ranj=max(0,v1):max.factorial(x[t-1],runif(1)) x[t]=sample(ranj,size=1) } barplot(as.vector(rbind( table(x)/length(x),dpois(min(x):max(x), lambda))),col=c("sienna","gold"))

**A**s you can easily check by running the code, it does not work. My student actually majored my MCMC class and he spent quite a while pondering why the code was not working. I did ponder as well for a part of a morning in Warwick, removing causes for exponential or factorial overflows (hence the shape of the code), but not eliciting the issue… (This now sounds like lethal fugu sashimi! ) Before reading any further, can you spot the problem?!

**T**he corrected R code is as follows:

x = rep(lambda, nsim) for (t in 2:nsim){ v1=ceiling((log(runif(1))/log(lambda))+x[t-1]) ranj=max(0,v1):max.factorial(x[t-1],runif(1)) if (length(ranj)>1){ x[t] = sample(ranj, size = 1) }else{ x[t]=ranj} }

**T**he culprit is thus the R function * sample* which simply does not understand Dirac masses and the basics of probability! When running

> sample(150:150,1) [1] 23

you can clearly see where the problem stands…! Well-documented issue with * sample* that already caused me woes… Another interesting thing about this slice sampler is that it is awfully slow in exploring the tails. And to converge to the centre from the tails. This is not very pronounced in the above graph with a mean of 6. Moving to 50 makes it more apparent:

**T**his is due to the poor mixing of the chain, as shown by the raw sequence below, which strives to achieve a single cycle out of 10⁵ iterations! In any case, thanks to Kévin for an interesting morning!

## sliced Poisson

Posted in Books, Kids, pictures, R, Running, Statistics, University life with tags ENSAE, Monte Carlo Statistical Methods, Poisson distribution, R, slice sampling on March 18, 2014 by xi'an**O**ne of my students complained that his slice sampler of a Poisson distribution was not working when following the instructions in Monte Carlo Statistical Methods (Exercise 8.5). This puzzled me during my early morning run and I checked on my way back, even before attacking the fresh baguette I had brought from the bakery… The following R code is the check. And it does work! As the comparison above shows…

slice=function(el,u){ #generate uniform over finite integer set mode=floor(lambda) sli=mode x=mode+1 while (dpois(x,el)>u){ sli=c(sli,x);x=x+1} x=mode-1 while (dpois(x,el)>u){ sli=c(sli,x);x=x-1} return(sample(sli,1))} #example T=10^4 lambda=2.414 x=rep(floor(lambda),T) for (t in 2:T) x[t]=slice(lambda,runif(1)*dpois(x[t-1],lambda)) barplot(as.vector(rbind( table(x)/length(x),dpois(0:max(x), lambda))),col=c("sienna","gold"))