# Statisfaction

## the much smaller world of pseudo-random generation

Posted in General by nicolaschopin on 16 September 2013

I start a new course at ENSAE this year, on “Monte Carlo and simulation methods”. I intend to cover pseudo-random generators at the beginning, so I’m thinking about how to teach this material which I’m not so familiar with.

One very naive remark: in a “truly random world”, when I flip a coin $n$ times, I obtain one out of $2^n$ possible outcomes, with probability $2^{-n}$. In the real world, if I use a computer to toss $n$ coins, the number of possible outcomes (for these $n$ successive tosses) is bounded by $2^{32}$. This is because a stream of pseudo-random numbers is completely determined by the seed (the starting point of the stream), and most generators are based on 32-bits seeds.

Compare $2^{32}$ with $2^n$ when $n$ is large, and you see that PRNG is quite a crude approximation of randomness. Of course, it’s not so bad in practice, because usually you are not interested in the exact value of a vector of $n$ successive coin tosses, but rather at some summary of dimension $d\ll 2^{32}$. Still, the pseudo-random world is much smaller than the random world it is supposed to mimic.

I found this remark quite scary, and I think I’ll use it to impress on my students the limitations of PRNG. By the way, if you like horror stories on PRNG, you might find  the slides of Régis Lebrun (for a talk at BigMC he gave a few years back) quite entertaining. It was really funny to see the faces of my colleagues turning white as Régis was giving more and more evidence that we are often too confident in PRN generators and oblivious of their limitations. I suspect my own face was very much the same colour.

Tagged with: ,

## A newcomer at Statisfaction

Posted in Statistics by nicolaschopin on 7 September 2013

I am Nicolas Chopin,  a Professor of Statistics at the ENSAE, and my colleagues and good friends that manage Statisfaction kindly agreed that I would join their blog. I work mostly on “Bayesian Computation”, i.e. Monte Carlo and non-Monte Carlo methods to compute Bayesian quantities; a strong focus of my research is on Sequential Monte Carlo (aka particle filters).

I don’t plan to blog very regularly, and only on stuff related to my research, at least in some way. Well, that’s the idea for now. Stay tuned!

Nicolas

## From SVG to probability distributions [with R package]

Posted in R, Statistics by Pierre Jacob on 25 August 2013

Hey,

To illustrate generally complex probability density functions on continuous spaces, researchers always use the same examples, for instance mixtures of Gaussian distributions or a banana shaped distribution defined on $\mathbb{R}^2$ with density function:

$f(x,y) = \exp\left(-\frac{x^2}{200} - \frac{1}{2}(y+Bx^2-100B)^2\right)$

If we draw a sample from this distribution using MCMC we obtain a [scatter]plot like this one:

Fig. 1: a sample from the very lame banana shaped distribution

Clearly it doesn’t really look like a banana, even if you use yellow to colour the dots like here. Actually it looks more like a boomerang, if anything. I was worried about this for a while, until I came up with a more realistic banana shaped distribution:

Fig. 2: a sample from the realistic banana shaped distribution

See how the shape is well defined compared to the first figure? And there’s even the little tail, that proves so convenient when we want to peel off the fruit. More generally we might want to create target density functions based on general shapes. For this you can now try RShapeTarget, which you can install directly from R using devtools:

library(devtools)


The package parses SVG files representing shapes, and creates target densities from them. More precisely, a SVG files contains “paths”, which are sequence of points (for instance the above banana is a single closed path). The associated log density at any point $x$ is defined by $-1/(2\lambda) \times d(x, P)$ where $P$ is the closest path of the shape from $x$ and $d(x,P)$ is the distance between the point and the path. The parameter $\lambda$ specifies the rate at which the density decays when the point goes away from the shape. With this you can define the maple leaf distribution, as a tribute to JSM 2013:

Fig. 3: a sample the “O Canada” probability distribution.

In the package you can get a distribution from a SVG file using the following code:

library(RShapeTarget)
# create target from file
my_shape_target <- create_target_from_shape(my_svg_file_name, lambda =1)
# test the log density function on 25 randomly generated points
my_shape_target$logd(matrix(rnorm(50), ncol = 2), my_shape_target$algo_parameters)


Since characters are just a bunch of paths, you can also define distributions based on words, for instance:

Hodor: Hodor.

which is done as follows (warning you’re only allowed a-z and A-Z, no numbers no space no punctuation for now):

library(RShapeTarget)
word_target <- create_target_from_word("Hodor")


For the words, I defined the target density function as before, except that it’s constant on the letters: so if a point is outside a letter its density is computed based on the distance to the nearest path; if it’s inside a letter it’s just constant, so that the letters are “filled” with some constant density. I thought it’d look better.

Now I’m not worried about the banana shaped distribution any more, but by the fact that the only word I could think of was “Hodor” (with whom you can chat over there).

## Joint Statistical Meeting 2013

Posted in General, Seminar/Conference, Statistics by Pierre Jacob on 23 July 2013

A typical statistical meeting.

Hey,

In a few weeks (August 3-8) I’ll attend the Joint Statistical Meeting in Montréal, Canada. According to Wikipedia it’s been held every year since 1840 and now gathers more than 5,000 participants!

I’ll talk in a session organized by Scott Schmidler, entitled Adaptive Monte Carlo Methods for Bayesian Computation; you can find the session programme here [online program]. I’ll talk about score and Fisher observation matrix estimation in state-space models.

According to the rumour and Christian’s reflections on the past years (2009, 2010, 2011), I should prepare my schedule in advance to really enjoy this giant meeting. So if you want to meet there, please send me an e-mail!

See you in Montréal!

## Path storage in the particle filter

Posted in Statistics by Pierre Jacob on 12 July 2013

Typical ancestry tree generated by a particle filter

Hey particle lovers,

With Lawrence Murray and Sylvain Rubenthaler we looked at how to store the paths in the particle filter, and the related expected memory cost. We just arXived a technical report about it. Would you like to know more?

## Intractable likelihoods, unbiased estimators and sign problem

Posted in Statistics by Pierre Jacob on 1 July 2013

Computing a Taylor expansion with random truncation can be done Swiftly.

Hey all,

We’re at the Big Data era blablabla, but the advanced computational methods usually don’t scale well enough to match the increasing sizes of datasets. For instance, even in a simple case of i.i.d. data $y_1, y_2, \ldots y_n$ and an associated likelihood function $\mathcal{L}(\theta; y_1, y_2, \ldots, y_n)$, the cost of evaluating the likelihood function at any parameter $\theta$ is typically growing at least linearly with $n$. If you then plug that likelihood into an optimization technique to find the Maximum Likelihood Estimate, or into a sampling technique such as Metropolis-Hastings to sample from the posterior distribution, the computational cost grows accordingly for a fixed number of iterations. However you can get unbiased estimates of the log-likelihood by drawing $m < n$ points $i_1, \ldots, i_m$ uniformly in the index set $\{1, \ldots, n\}$ and by computing $(n/m) \log \mathcal{L}(\theta; y_{i_1}, \ldots, y_{i_m})$. This way you sub-sample from the whole dataset, and you can choose $m$ according to your computational budget. However is it possible to perform inference with these estimates instead of the complete log-likelihood?

## Call for contributors

Posted in General by Pierre Jacob on 24 June 2013

Hey,

This blog started as a collaborative blog written by then PhD students at CREST. Now some of us have left the lab but still feel like blogging from time to time so we use this blog. Going further, I don’t see any reason not to broaden our perspective by letting other people participate, in order to maintain a decent activity. The target would be at least a post per week. I am sure that many junior researchers out there feel like they could write a post or two, so this could be the place to share your views!

If you’re interested, either in a one-time blog post or on a more regular basis, please feel free to contact us by e-mail or in the comments below. You can just browse through the blog if you’re not sure what the scope is… actually the scope is pretty ill-defined, but includes tips and tricks in R, LaTeX, conferences in Statistics, mostly Bayesian or computational, random datasets, recent articles, reports on unusual use of statistical methods…

Cheers!

## PhDTree and juvenile fishes

Posted in General by Pierre Jacob on 5 June 2013

Hey,

I’ve just discovered this website called PhDTree, which is awesome because it shows that science-wise I am a great-great-great-great-great-great-great child of Siméon Denis Poisson.

And then I’ve looked at the genealogy of all of my French colleagues and they are all descendants of him, which is rather annoying.

## Derivative-free estimate of derivatives

Posted in Statistics by Pierre Jacob on 23 April 2013

Two Hessian soldiers having a ball.

Hey,

Arnaud Doucet, Sylvain Rubenthaler and I have just put a technical report on arXiv about estimating the first- and second-order derivatives of the log-likelihood (also called the score and the observed information matrix respectively) in general (intractable) statistical models, and in particular in (non-linear non-Gaussian) state-space models. We call them “derivative-free” estimates because they can be computed even if the user cannot compute any kind of derivatives related to the model (as opposed to e.g. this paper and this paper). Actually in some cases of interest we cannot even evaluate the log-likelihood point-wise (we do not have a formula for it), so forget about explicit derivatives. Would you like to know more?

## BayesComp on wikidot

Posted in General by Pierre Jacob on 16 April 2013

Since the classic picture of Bayes is actually not a picture of Bayes, we might as well use Ryan.

Hey,

Just a quick note about BayesComp, a new wiki about Bayesian Computational Statistics (see this outdated but well-written introduction if you really don’t know what that is), as Xian pointed out. It is organised by the ISBA Section on Bayesian Computation, notably Peter Green and Nicolas Chopin so far. If the community gets into it, it could become the nerve centre for online resources about Bayesian Computation, which so far are quite scattered and poorly advertised.

Good luck to BayesComp!