## fiducial simulation

Posted in Books, Kids, pictures, Statistics, Travel, University life with tags , , , , , , , , , , , on April 19, 2018 by xi'an

While reading Confidence, Likelihood, Probability), by Tore Schweder and Nils Hjort, in the train from Oxford to Warwick, I came upon this unexpected property shown by Lindqvist and Taraldsen (Biometrika, 2005) that to simulate a sample y conditional on the realisation of a sufficient statistic, T(y)=t⁰, it is sufficient (!!!) to simulate the components of  y as y=G(u,θ), with u a random variable with fixed distribution, e.g., a U(0,1), and to solve in θ the fixed point equation T(y)=t⁰. Assuming there exists a single solution. Brilliant (like an aurora borealis)! To borrow a simple example from the authors, take an exponential sample to be simulated given the sum statistics. As it is well-known, the conditional distribution is then a (rescaled) Beta and the proposed algorithm ends up being a standard Beta generator. For the method to work in general, T(y) must factorise through a function of the u’s, a so-called pivotal condition which brings us back to my post title. If this condition does not hold, the authors once again brilliantly introduce a pseudo-prior distribution on the parameter θ to make it independent from the u’s conditional on T(y)=t⁰. And discuss the choice of the Jeffreys prior as optimal in this setting even when this prior is improper. While the setting is necessarily one of exponential families and of sufficient conditioning statistics, I find it amazing that this property is not more well-known [at least by me!]. And wonder if there is an equivalent outside exponential families, for instance for simulating a t sample conditional on the average of this sample.

## Le Monde puzzle [#1049]

Posted in Books, Kids, R with tags , , , on April 18, 2018 by xi'an

An algorithmic Le Monde mathematical puzzle with a direct

Alice and Bob play a game by picking alternatively one of the remaining digits between 1 and 10 and putting it in either one of two available stacks, 1 or 2. Their respective gains are the products of the piles (1 for Alice and 2 for Bob).

The problem is manageable by a recursive function

facten=factorial(10)
pick=function(play=1,remz=matrix(0,2,5)){
if ((min(remz[1,])>0)||(min(remz[2,])>0)){#finale
remz[remz==0]=(1:10)[!(1:10)%in%remz]
return(prod(remz[play,]))
}else{
gainz=0
for (i in (1:10)[!(1:10)%in%remz]){
propz=rbind(c(remz[1,remz[1,]>0],i,
rep(0,sum(remz[1,]==0)-1)),remz[2,])
gainz=max(gainz,facten/pick(3-play,remz=propz))}
for (i in (1:10)[!(1:10)%in%remz]){
propz=rbind(remz[1,],c(remz[2,remz[2,]>0],i,
rep(0,sum(remz[2,]==0)-1)))
gainz=max(gainz,facten/pick(3-play,remz=propz))}
return(gainz)}}


that shows the optimal gain for Alice is 3360=2x5x6x7x 8, versus Bob getting 1080=1x3x4x9x10. The moves ensuring the gain are 2-10-…

Posted in Kids, pictures, Travel with tags , , , , , , , on April 14, 2018 by xi'an

I found out most recently that it was conceivable to make one own’s soba noodles at home! (Soba means buckwheat.) And hence embarked on the adventure on two consecutive weekends (one batch of noodles each!). As I like very much these particular noodles, having trouble finding manufactured ones around (the last ones I bought had no buckwheat flour content and were made in Hungary…). The recipe is awfully simple, since it consists in making a dough out of buckwheat and wheat flours, flattening it out (and again), and cutting thin noodles from the resulting folds. In practice, it does not work so well, from the dough crumbling when getting very thin (as buckwheat does not have the same adherence as wheat, missing the gluten component), to the folds sticking with one another when preparing to cut four or eight layers at a time. So far the results are more tagliatelle like than soba, but still eatable as is (especially for lunch at the office) and giving me a challenge to try next new techniques (like using hot water when making the dough).

## Masterclass in Bayesian Statistics in Marseilles next Fall

Posted in Books, Kids, Mountains, pictures, R, Running, Statistics, Travel, University life with tags , , , , , , , , , , , , , , , , on April 9, 2018 by xi'an

This post is to announce a second occurrence of the exciting “masterclass in Bayesian Statistics” that we organised in 2016, near Marseilles. It will take place on 22-26 October 2018 once more at CIRM (Centre International de Recherches Mathématiques, Luminy, Marseilles, France). The targeted audience includes all scientists interested in learning how Bayesian inference may be used to tackle the practical problems they face in their own research. In particular PhD students and post-docs should benefit most directly from this masterclass. Among the invited speakers, Kerrie Mengersen from QUT, Brisbane, visiting Marseilles this Fall, will deliver a series of lectures on the interface between Bayesian statistics and applied modelling, Havard Rue from KAUST will talk on computing with INLA, and Aki Vehtari from Aalto U, Helsinki, will give a course on Bayesian model assessment and model choice. There will be two tutorials on R and on Stan.

All interested participants in this masterclass should pre-register as early as possible, given that the total attendance is limited to roughly 90 participants. Some specific funding for local expenses (i.e., food + accommodation on-siteat CIRM) is available (thanks to CIRM, and potentially to Fondation Jacques Hadamard, to be confirmed); this funding will be attributed by the scientific committee, with high priority to PhD students and post-docs.

## truncated Gumbels

Posted in Books, Kids, pictures, Statistics with tags , , , , , , , on April 6, 2018 by xi'an

As I had to wake up pretty early on Easter morning to give my daughter a ride, while waiting I came upon this calculus question on X validated of computing the conditional expectation of a Gumbel variate, conditional on its drifted version being larger than another independent Gumbel variate with the same location-scale parameters. (Just reminding readers that a Gumbel G(0,1) variate is a double log-uniform, i.e., can be generated as X=-log(-log(U)).) And found after a few minutes (and a call to Wolfram Alpha integrator) that

$\mathbb{E}[\epsilon_1|\epsilon_1+c>\epsilon_0]=\gamma+\log(1+e^{-c})$

which is simple enough to make me wonder if there is a simpler derivation than the call to the exponential integral Ei(x) function. (And easy to check by simulation.)

Incidentally, I discovered that Emil Gumbel had applied statistical analysis to the study of four years of political murders in the Weimar Republic, demonstrating the huge bias of the local justice towards right-wing murders. When he signed the urgent call [for the union of the socialist and communist parties] against fascism in 1932, he got expelled from his professor position in Heidelberg and emigrated to France, which he had to leave again for the USA on the Nazi invasion in 1940. Where he became a professor at Columbia.

## Springer no more!

Posted in Books, Kids, Statistics, University life with tags , , , , , on April 4, 2018 by xi'an

Just learned that, starting from tomorrow night, I will not have access to any of the Springer journals, as the negotiations between the consortium of French universities, research institutes, higher educations schools, and museums, failed. The commercial published refusing to stem the ever increasing fees, while happily taking in the fast increasing open access fees it pressures from authors, a unique example of triple taxation (researchers’ salaries, open access duties, and enormous non-negotiable subscription rates for the whole package of journals)… Following their German counterparts. Well, this is an opportunity for the boards of all these journals to withdraw and create the phantom version of their formal journal, evaluating and reviewing papers already available on arXiv! And I should definitely get my acts together, rise from my winter-is-coming lethargy, and launch PCI Comput Stat now!!!

## Le Monde puzzle [#1048]

Posted in Books, Kids, R with tags , , , , , on April 1, 2018 by xi'an

An arithmetic Le Monde mathematical puzzle:

A magical integer m is such that the remainder of the division of any prime number p by m is either a prime number or 1. What is the unique magical integer between 25 and 100? And is there any less than 25?

The question is dead easy to code

primz=c(1,generate_primes(2,1e6))
for (y in 25:10000)
if (min((primz[primz>y]%%y)%in%primz)==1) print(y)


and return m=30 as the only solution. Bon sang but of course!, since 30=2x3x5… (Actually, the result follows by dividing the quotient of the division of a prime number by 2 by 3 and then the resulting quotient by 5: all possible cases produce a remainder that is a prime number.) For the second question, the same code returns 2,3,4,6,8,12,18,24 as further solutions. There is no solution beyond 30.