particles 0.3: waste-free SMC, Fortran dependency removed, binary spaces

I have just released version 0.3 of particles (my Python Sequential Monte Carlo library). Here are the main changes:

No more Fortran dependency

Previous versions of particles relied on a bit of Fortran code to produce QMC (quasi-Monte Carlo) points. This code was automatically compiled during the installation. This was working fine for most users, but not all, unfortunately.

The latest (1.7) version of Scipy includes a stats.qmc sub-module. Particles 0.3 relies on this sub-module to generate QMC points, and thus is a pure Python package. This should mean fewer headaches when installing particles. Please let me know if this new version is indeed easier to install for you. Of course, make sure you have updated Scipy before installing particles; e.g. conda update scipy if you are using conda.

Waste-free SMC

With Dang, we wrote a paper on a new class of SMC samplers, waste-free SMC; see this paper on arxiv (to be published soon in JRSSB). In particular, the new version describes a particular scenario where it is possible to show formally that waste-free SMC >> standard SMC (in the sense of lower asymptotic variance).

The module smc_samplers now implements waste-free SMC by default (but standard SMC is still available, through option wastefree=False). Check the following notebook to see how to run an SMC sampler in particles.

SMC samplers for binary spaces

The new module binary_smc implements SMC samplers for binary spaces, i.e. {0, 1}^d, following Chopin and Schäfer (2014).

papers folder

The package now includes a folder called “papers”, which contains scripts that reproduce selected numerical experiments from previous papers:

  • scripts in sub-folder binarySMC reproduce most of the numerical experiments from Schäfer and Chopin (2014).
  • a script in sub-folder wastefreeSMC reproduces the first numerical experiment of Dau & Chopin (2020) on logistic regression. (See Dang’s github repo for the other experiments.)


  • Added a new resampling scheme, called killing (which may be traced to papers and work by Pierre del Moral).
  • Added a tutorial notebook on how to define non-trivial state-space models.

Comments / questions / poems ?

If you want to try particles, the first thing to read is the notebook tutorials. Second thing is to read the documentation of the respective modules. If you are still lost, feel free to raise an issue on github (or send me an e-mail, but github issues are more practical).

Online workshop: Measuring the quality of MCMC output

Yes, we need a better logo.

Hi all,

With Leah South from QUT we are organizing an online workshop on the topic of “Measuring the quality of MCMC output”. The event website is here with more info:

This is part of ISBA BayesComp section’s efforts to organize activities while waiting for the next “big” in-person meeting, hopefully in 2023. The event benefits from the generous support of QUT Centre for Data Science. The event’s website will be regularly updated between now and the event in October 2021, with three live sessions:

  • 11am-2pm UTC on Wednesday 6th October,
  • 1pm-4pm UTC on Thursday 14th October,
  • 3pm-6pm UTC on Friday 22nd October.

Registration is free but compulsory (form here) as we want to make sure the live sessions remain convivial and focused; hence the rather specific theme, but it’s an exciting topic with lots of very much open questions, which we hope will attract both practitioners and methodologists. Meanwhile some material will be available on the website to everyone, including video recordings of presentations, and posters, so that the workshop hopefully benefits the wider community.

If you have suggestions for this event, or would like to organize a similar event in the future, on another “BayesComp” topic, do not hesitate to get in touch. Our contact details are on the workshop’s website.

Dempster’s analysis and donkeys

This post is about estimating the parameter of a Bernoulli distribution from observations, in the “Dempster” or “DempsterShafer” way, which is a generalization of Bayesian inference. I’ll recall what this approach is about, and describe a Gibbs sampler to perform the computation. Intriguingly the associated Markov chain happens to be equivalent to the so-called “donkey walk” (not this one), as pointed out by Guanyang Wang and Persi Diaconis.

Continue reading “Dempster’s analysis and donkeys”

particles 0.2: what’s new, what’s next (your comments most welcome)

I have just released version 0.2 of my SMC python library, particles. I list below the main changes, and discuss some ideas for the future of the library.

New module: variance_estimators

This module implements various variance estimators that may be computed from a single run of an SMC algorithm, à la Chan and Lai (2013) and Lee and Whiteley (2018). For more details, see this notebook.

New module: datasets

This module makes it easier to load the datasets included in the module. Here is a quick example:

from particles import datasets as dts

dataset = dts.Pima()
help(dataset) # basic info on dataset
help(dataset.preprocess) # how data was pre-processed
data = # typically a numpy array

Performance improvements: multiprocessing

The library makes it possible to run several SMC algorithms in parallel, using the multiprocessing module. Hai-Dang Dau noticed there was some performance issue with the previous implementation (a few cores could stay idle) and fixed it.

While testing the new version, I noticed that function distinct_seeds (module utils), which, as the name suggests, generate distinct random seeds for the processes run in parallel, could be very slow in certain cases. I changed the way the seeds were generated to fix the problem (using stratified resampling). I will discuss this in more detail a separate blog post.

Plans for the future

Development of this library is partly driven by interactions with users. For instance, the next version will have a more general MvNormal distribution (allowing for a covariance matrix that varies across particles), because one colleague got in touch and needed that feature.

So don’t be shy, if you don’t see how to do something with particles, please get in touch. It’s likely our interaction will help me to either improve the documentation or add new, useful features. Of course, I also welcome direct contributions (through pull requests)!

Otherwise, I have several ideas for future releases, but, for the next one, it is likely I will focus on the following two areas.

SMC samplers

My priority #1 is to implement waste-free SMC in the package, following our recent paper with Dang. (Dang already has released his own implementation, which is built on top of particles, but, given that waste-free SMC seems to offer better performance than standard SMC samplers, it seems important to have it available in particles).

When this is done, I plan to add several important applications of SMC samplers, such as:

I also plan to document SMC samplers a bit better.

integration with Pytorch (or JAX, Tensorflow, etc)

Python libraries such as Tensorflow, Pytorch or JAX are all the rage in machine learning. They offer access to very fancy stuff, such as auto-differentation, and computation on the GPU.

I have started to play a bit with Pytorch, and even have a working implementation of a particle filter that runs entirely on the GPU. The idea is to make the core parts of particles completely independent of numpy. In that way, one may use Pytorch tensors to store the particles and their weights. This is really work in progress.

On the benefits of reviewing papers

Would you have agreed to review this paper if you had been asked?

When I’m asked by students whether they should accept some referee invitation (being it for a stat journal or a machine learning conference) I almost invariably say yes. I think that there is a lot to be learnt when refereeing papers and that this worth the time spent in the process. I’ll detail in this post why I think so.

First, this post is not about tips on how to write a referee report, but rather on why. It is instructive to consult tips on the hows, and good posts can be found out there. Note that some journals will also have specific guidelines.

Before diving into the benefits of refereeing, let me first say that a referee invitation can also be declined for many good reasons: in case of a conflict of interest (CoI), and/or if some of the authors are too close to you in some sense (although in some fields with a tiny community, this almost inevitably happens); if you do not feel qualified enough; or sometimes, if you feel you are qualified, but the refereeing task can seem overwhelming due to length or technicality of the paper; do not feel obliged to accept invitations from journals you do not know about, and of course ignore those coming from predatory journals or publishers (use this checklist). In any case, be conscious that it is ok to decline an invitation. Keep in mind that the associate editor in charge will very much appreciate pointers to alternative referee names.

Now, what are the benefits of refereeing? It is a legitimate question, given that refereeing work is usually time-consuming, done on a voluntary basis, without implying any direct or instant reward. So it is important to understand what you can gain out of it.

Learning about editorial process

In the early stages of an academic career, refereeing papers is an opportunity of learning by doing about the ins and outs of the editorial mechanism. You do not get the chance to practice replying to referee reports every other day when you are a student. But getting papers to review, you may also get to see replies by the authors, and reports from other referees (eg in revision rounds). This may help and build some habit about how you will get into action when your turn comes to reply to referees!

Opening research interests

We are usually asked to referee papers in our own area of expertise, but accepting to review papers slightly outside of one’s research interests can be rewarding. Be curious! There is a chance that reading submitted papers will trigger new research directions of yours. This happened to me at least twice: I have started to work on Bayesian deep learning after refereeing an ICLR paper dealing with the behaviour of neural networks in the infinitely wide limit; and (dis)proving a conjecture stated in a COLT submission stimulated a new line of research of mine on the sub-Gaussian property of random variables. Pay attention that in order to start working on such submitted papers in a legit way, you should ensure that they are also made available as preprints on some open repository like arxiv.

Prompting new opportunities

Refereeing papers surely increases your visibility. It is also a preliminary step before being associate editor. I’m AE for several stat journals, and managing papers is a task that I find enjoyable, with a social side that consists of writing referee invitation messages to colleagues. This helps connect or stay in touch with colleagues we do not have occasions to meet in conferences those days!

Everything You Always Wanted to Know About SMC, but were afraid to ask

Ever wanted to learn more about particle filters, sequential Monte Carlo, state-space/hidden Markov models, PMCMC (particle MCMC) , SMC samplers, and related topics?

In that case, you might want to check the following book from Omiros Papaspiliopoulos and I, which has just been released by Springer:

and which may be ordered from their web-site, or from your favourite book store.

The aim of the book is to cover the many facets of SMC: the algorithms, their practical uses in different areas, the underlying theory, how they may be implemented in practice, etc. Each chapter contains a “Python corner” which discusses the practical implementation of the covered methods in Python, a set of exercises, and bibliographical notes. Speaking of chapters, here is the table of contents:

  1. Introduction
  2. Introduction to state-space models
  3. Beyond state-space models
  4. Introduction to Markov processes
  5. Feynman-Kac models: definition, properties and recursions
  6. Finite state-spaces and hidden Markov models
  7. Linear-Gaussian state-space models
  8. Importance sampling
  9. Importance resampling
  10. Particle filtering
  11. Convergence and stability of particle filters
  12. Particle smoothing
  13. Sequential quasi-Monte Carlo
  14. Maximum likelihood estimation of state-space models
  15. Markov chain Monte Carlo
  16. Bayesian estimation of state-space models and particle MCMC
  17. SMC samplers
  18. SMC^2, sequential inference in state-space models
  19. Advanced topics and open problems

And here is one fancy plot taken from the book. (For some explanation, you will have to read it!)

A big thanks to all the colleagues who took the time to read draft versions and send feedback (see the introduction for a list of names). Also, don’t write books, folks. Seriously, it takes WAY too much time…

Categorical distribution, structure of the second kind and Gumbel-max trick

Hi all,

This post is about a way of sampling from a Categorical distribution, which appears in Arthur Dempter‘s approach to inference as a generalization of Bayesian inference (see Figure 1 in “A Generalization of Bayesian Inference”, 1968), under the name “structure of the second kind”. It’s the starting point of my on-going work with Ruobin Gong and Paul Edlefsen, which I’ll write about on another day. This sampling mechanism turns out to be strictly equivalent to the “Gumbel-max” trick that got some attention in machine learning see e.g. this blog post by Francis Bach.

Continue reading “Categorical distribution, structure of the second kind and Gumbel-max trick”

Course on Couplings and Monte Carlo

Hi everyone,

This short post is just to point to a course on “Couplings and Monte Carlo”, available here Versions of the course were given in Université Paris-Dauphine in February 2020 (thanks Robin Ryder and Christian P. Robert), at the University of Bristol in March 2020 (thanks Anthony Lee) and at the University of Torino for the M.Sc. in Stochastics and Data Science in May 2020 (thanks Matteo Ruggiero). I am grateful to these colleagues and their institutions for supporting this course. The course website points to about 100 pages of lecture notes, and 16 videos are available on youtube. It is intended for advanced undergraduate students or graduate students, with some previous exposure to Monte Carlo methods. This is work in progress, and as I am hoping to develop the course over the coming years, feedback would be much welcome.

Follow-up on my previous post on covid deaths in France

In my previous post, I compared two sources of data regarding death counts (one from SPF, Santé Publique France, for covid deaths in hospitals; one from INSEE, for all-cause deaths), in order to get a better idea of the actual death toll of covid-19 in France.

In this post, I would like to do the following:

  1. describe briefly a new, richer data-set recently published by INSEE (and do some graphs);
  2. use the updated data (from both sources) to repeat my analysis, with some variants (weekly aggregates, separating men and women);
  3. reply to a few comments I got on LinkedIn and elsewhere;
  4. provide a few pointers regarding death counts in other countries (particularly the UK).

New INSEE data

INSEE now provides every Friday an exhaustive data-set that records, for each death that has occurred since 01-01-2018, the following variables: date of birth, date of death, sex, département of death, and so on. Neat. Let’s take this opportunity to do a few plots, such as this one:

Number of deaths in France each day in 2018, 2019, 2020, until 13th Apr 2020 (as recored in INSEE dataset described above)

(it’s nice to observe this sharp drop) or that one:

Death counts during weeks 13, 14 and 15 per age, for years 2018-20 (same INSEE dataset

The latter plot covers the same period (weeks 13 to 15, 23rd March to 12th April) as in the analysis below. As expected, over-mortality seems to affect mostly people above 60.

Updated analysis

Ok, now let’s repeat my previous analysis, based on merging the SPF data (daily covid death counts in hospitals, in each département and each sex) and the aforementioned INSEE data (all-cause deaths). Except this time:

  1. The overlap between the two datasets now covers more than three weeks (18th March, first date in SPF dataset, to 12th April, latest date in INSEE dataset) so I decided to consider weekly aggregates, for two reasons: they are more stable than daily aggregates, and less affected by artifacts such as delays (e.g. a death occurring during a week-end is reported to the next Monday).
  2. I also separated men and women.
  3. I am going to simplify a bit the model, and simply regress excess deaths (number of deaths in 2020 minus the average over 2018 and 2019) on hospital deaths.

First, a joint plot:

So, to recap, each point in this plot corresponds to a pair of death counts, for each département in France, each week between 13 and 15, and for each sex. The corresponding linear regression (without an intercept) gives a slope estimate of 1.79 (95% confidence interval: [1.73, 1.85]). The basic interpretation would be: in each département, when 100 covid deaths occur in hospitals, the number of covid-related (see below) deaths should be approximately 179. The current total number of covid deaths reported by SPF is 22, 614, which is 60% above that the number of covid deaths in hospitals (14050). So this estimate suggests the actual death toll might be a tad larger. More about the interpretation below.

Now for something more interesting: let’s redo the previous plot, but with a different colour for each sex:

Clearly the two linear trends are different; see below the OLS estimates.

sexslope estimateslope 95% confidence intervalR^2
F2.40[2.30, 2.50]89%
M1.56[1.50, 1.62]90%

What is going on? Well, women tend to live longer than men. And the proportion of women in EHPADS (French retirement homes) is 74%. Since the main reason behind the discrepancy between hospital deaths and excess deaths is covid death occurring in pension homes, these results make sense.

Reply to previous comments

What’s the point?

Fair enough, since 4th April, SPF does include in its total estimate both hospital deaths and pension home deaths, and the proportion of the latter is not too far from my estimate. Note that however that:

  • it’s really hard to estimate properly the number of covid deaths occurring in pension homes. Apparently several pension homes did not provide any data, while others marked as “covid” all the deaths that have occurred after the first covid deaths.
  • My estimate might measure other direct or indirect effects of the pandemic, such as people dying at home, people not receiving proper care because the health system is at capacity and so on.
  • The fact that data from two different institutions may be compared, and seem to be somehow consistent, is, in my opinion, a good piece of news which deserves to be reported!

Accounting for car accidents

Boy, that one was popular. Please have a look a the plot on the front page of ONIRS (click on “tués”): yes, the number car-related deaths dropped sharply thanks to the lock-down… But in March of last year, this number was around 250, that is, 1% of the current covid death count. “Fun” fact: this point would have been quite relevant in the 70s! In those years, the number of car-related deaths was about five times larger (18 034 deaths in 1972).

Better predictive model

The idea of comparing the 2020 deaths to the average of the two previous years is a bit crude, and demographers have better models to predict death counts based on age repartition and so on. That said, the notion of “excess deaths” seems quite popular in various countries, as I explain below, so I guess that my approach is not so daft after all.

In other countries

To be honest, I was hoping to apply the same approach to the UK, a country where the official estimate is still limited to hospital deaths, and thus clearly quite biased; see e.g. this Guardian paper. Sadly, Public Health England only reports daily hospital death counts … per nation (nation=England, Scotland, Wales, or Northern Ireland). On the other hand, the Office of National Statistics reports every week the number of “excess deaths” (relative to the five year average), and the proportion of these deaths where the word “covid” is mentioned on the death certificate.

Interestingly, the Guardian paper I mentioned above first complains that the UK only reports hospital deaths, and then claims erroneously that the UK is still behind France in terms of covid mortality. It’s not, if you compare in terms of hospital deaths (UK: 20,319 on Saturday, while France: 14,050). The fact that even journalists reporting on this issue may get it wrong seems indicative of how confusing are covid death data.

More generally, my impression is that looking at “excess deaths” makes far more sense for most countries at the moment: it’s easier to measure (albeit with a delay of course), and easier to interpret. This is also more or less the point made by this NYT paper. (Notice how their plot for France only covers January to April; for the complete plot, see my first plot above!).

%d bloggers like this: