Le Monde puzzle [52]

The last puzzle of the year in Le Monde reads as follows (as far as I understand its wording!):

Iter(n,x,y) is the function

Iter=function(n,x,y){

 if (n==1){
 output=trunc(y/10)+x*(y%%10)
 }else{
 output=Iter(n-1,x,Iter(1,x,y))}

 return output
 }

Find the seven-digit number z such that
Iter(6,1,z)=12, Iter(6,2,z)=19, Iter(6,3,z)=29,
and Iter(6,-1,z)=Iter(6,-2,z)=Iter(6,-3,z)=0.

Obviously, the brute-force solution of listing all 90 million seven digit numbers until the six constraints are met is feasible (especially around New Year since the mainframe computer is completely at rest!). However, this sounds like the last resort solution and I thus tried first a simulated annealing approach already tested for the sudoku problem a few years ago… (This puzzle is actually of the same nature as the sudoku problem,  in particular because we do know when we find the solution, except that checking for the six conditions to hold is apparently not so straightforward. For us if not for the computer.)

the half-made world

“He was–ah–he was carving the roast at dinner, for the Dean of the Faculty of Mathematics and the Bishop of Lodenstein, and others. And the effort, and the occasion, were too much for him. He swelled up in his shirtsleeves and burst. He fell with his mustaches in the gravy.”

tor.com put the first chapter of this book on-line and it was original enough to make me buy it… I read the remainder of the half-made world over the past week and really enjoyed the weirdness of it! It came as part of the steampunk month on tor.com but it blurs the genre boundaries as it also qualifies as weird West…

Posters at Adap[/MCM]’skiii

Due to a very small number of (declared) submissions—contrasting with the stable number of attendees—for the poster session on the evening of January 4 at Adap’skii, I have decided to cancel this session. The posters will be presented on the evenings of Wednesday, January 5 and Thursday, January 6, during MCMSki III. This will ensure an adequate audience and level of interaction for the poster presenters.

Christophe Andrieu being unfortunately unavailable at the present time, I am also unable to update the webpage prior to the meeting. Besides the above cancellation, the other changes in the program are

• Jan. 3, 6pm: the discussant of Faming Liang’s talk is Pierre Jacob (Université Paris Dauphine & CREST)
• Jan. 4, 8:30am: the discussant of Yves Atchadé is Luke Bornn (UBC, Vancouver)

Looking forward to the meetings in a few days! And to the 3 feet of fresh snow!!!

More typos in Chapter 5

Following Ashley’s latest comments on Chapter 5 of Introducing Monte Carlo Methods with R, I realised Example 5.5 was totally off-the-mark! Not only the representation of the likelihood should have used prod instead of mean, not only the constant should call the val argument of integrate, not only integrate  uses lower and upper rather than from and to, but also the approximation is missing a scale factor of 10, squared root of the sample size… The corrected R code is thus

> cau=rcauchy(10^2)
> mcau=median(cau)
> rcau=diff(quantile(cau,c(.25,.75)))/sqrt(10^2)
> f=function(x){
+   z=dcauchy(outer(x,cau,FUN="-"))
+   apply(z,1,prod)}
> fcst=integrate(f,lower=-20,upper=20)$val
> ft=function(x){f(x)/fcst}
> g=function(x){dt((x-mcau)/rcau,df=49)/rcau}
> curve(ft,from=-1,to=1,xlab="",ylab="",lwd=2)
> curve(g,add=T,lty=2,col="steelblue",lwd=2)

and the corrected Figure 5.5 is therefore as follows. Note that the fit by the t distribution is not as perfect as before. A normal approximation would do better.

This mistake is most embarrassing and I cannot fathom how I came with this unoperating program! (The more embarrassing as Cauchy‘s house is about 1k away from mine…) I am thus quite grateful to Ashley for her detailed study of this example.

Morning light

nlm [unused argument(s) (iter = 1)]

Ashley put the following comment on Chapter 5 of Introducing Monte Carlo Methods with R”:

I am reading chapter 5. I try to reproduced the result on page 128. The R codes don’t work on my laptop. When I try to run the following codes on page 128

> for (i in 1:(nlm(like,sta)$it)){
+ mmu=rbind(mmu,nlm(like,sta,iter=i)$est)}

I always get the error message

Error in f(x, …) : unused argument(s) (iter = 1)

It seems that the nlm function doesn’t accept the argument iter. I don’t know how to deal with it. I am in US. I guess the nlm version available to US R users is different from the version in EU. Please help.

And indeed with the most recent versions of R, like 2.12.1 on my own machine, calling nlm

> args(nlm)
function (f, p, ..., hessian = FALSE, typsize = rep(1, length(p)),
 fscale = 1, print.level = 0, ndigit = 12, gradtol = 1e-06,
 stepmax = max(1000 * sqrt(sum((p/typsize)^2)), 1000), steptol = 1e-06,
 iterlim = 100, check.analyticals = TRUE)

with the abbreviated argument iter instead of iterlim produces the above error message. This means the full syntax iterlim=i should now be used. In addition, the function nlm produces the minimum of the first argument f and like should thus be defined as

> like=function(mu){
+   -sum(log((.25*dnorm(da-mu[1])+.75*dnorm(da-mu[2]))))}

to end up with local maxima as on Figure 5.2. (Note: I do not think there are US versus EU versions of R…)

Ashley also pointed out another mistake on that page, namely that we used

> da=rbind(rnorm(10^2),2.5+rnorm(3*10^2))

instead of

> da=c(rnorm(10^2),2.5+rnorm(3*10^2))

to create a sample. Since the two normal samples have different sizes, rbind induces a duplication of the smaller sample, not what we intended!

https://xianblog.wordpress.com/2010/09/09/typo-in-chapter-5/comment-page-1/#comment-7677

Poster at MCMSki III

Here is the poster presented at MCMSki III next week by Pierre Jacob about our joint paper on parallelisation: