Archive for permutations

mostly Monte Carlo, the return²⁵

Posted in pictures, Statistics, University life with tags , , , , , , , , , , , , , , on October 9, 2025 by xi'an

Our local Mostly (and monthly) Monte Carlo seminar is back for a new academic year, now organized by Antoine Luciano and Timothy Johnston. The first session will take place at the PariSanté Campus on Friday 17 October 2025 (3:00pm, room 07), with the organisers opening the dance, with two talks:

3pm Timothy Johnston (CEREMADE, Université Paris Dauphine–PSL): Differential Privacy of Markov Chains

Joint work with Andrea Bertazzi, Alain Durmus and Gareth Roberts

In this talk we shall discuss differential privacy, a framework for quantifying the extent to which a random output depends on the information used to produce it. After introducing several related definition of differential privacy, we shall discuss techniques used to show the differential privacy of both trajectories and single draws from Markov Chains. In doing so we shall touch on a perturbation technique which allows for Wasserstein type bounds to be converted into stronger distances like the KL and Renyi divergence.

4pm Antoine Luciano (CEREMADE, Université Paris Dauphine–PSL): Permutations accelerate Approximate Bayesian Computation

Joint work with Charly Andral, Christian P. Robert and Robin J. Ryder

Approximate Bayesian Computation (ABC) methods have become essential tools for performing inference when likelihood functions are intractable or computationally prohibitive. However, their scalability remains a major challenge in hierarchical or high-dimensional models. In this paper, we introduce permABC, a new ABC framework designed for settings with both global and local parameters, where observations are grouped into exchangeable compartments. Building upon the Sequential Monte Carlo ABC (ABC-SMC) framework, permABC exploits the exchangeability of compartments through permutation-based matching, significantly improving computational efficiency. We then develop two further, complementary sequential strategies: Over Sampling, which facilitates early-stage acceptance by temporarily increasing the number of simulated compartments, and Under Matching, which relaxes the acceptance condition by matching only subsets of the data. These techniques allow for robust and scalable inference even in high-dimensional regimes. Through synthetic and real-world experiments – including a hierarchical Susceptible-Infectious-Recover model of the early COVID-19 epidemic across 94 French departments – we demonstrate the practical gains in accuracy and efficiency achieved by our approach.

permutations accelerate ABC!

Posted in Books, Kids, pictures, Statistics, University life with tags , , , , , , , , , , , , , , , , , , , , , , on July 9, 2025 by xi'an

Yesterday a arXival by Antoine Luciano, Charly Andral (both PhD students, now or then, at Paris Dauphine), Robin Ryder (formerly at Paris Dauphine, now at Imperial College London) and myself got posted. It proposes to improve the scalability of ABC methods by exploiting the (full or partial) exchangeability in the data by implementing permutation-based matching between observed and simulated samples. This significantly improves computational efficiency, which is further enhanced by sequential strategies such as over-sampling, which facilitates early-stage acceptance by temporarily increasing the number of simulated compartments, and under-matching, which relaxes the acceptance condition by matching only subsets of the data. The map of France appears in connection with an application of the method to estimating SIR parameters, department by department. (It is also reminding me of the cover of Markov Chain Monte Carlo methods in practice, the 1996 contributed book edited by Wally Gilks, Sylvia Richardson and David Spiegelhalter.)

overlap, overstreched

Posted in Books, Kids, R, Statistics with tags , , , , , , on June 15, 2020 by xi'an

An interesting challenge on The Riddler on the probability to see a random interval X’ing with all other random intervals when generating n intervals from Dirichlet D(1,1,1). As it happens the probability is always 2/3, whatever n>1, as shown by the R code below (where replicate cannot be replaced by rep!):

qro=function(n,T=1e3){
  quo=function(n){
     xyz=t(apply(matrix(runif(2*n),n),1,sort))
  sum(xyz[,1]<min(xyz[,2])&xyz[,2]>max(xyz[,1]))<0}
  mean(replicate(quo(n),T))}

and discussed more in details on X validated. As only a property on permutations and partitions. (The above picture is taken from this 2015 X validated post.)

Le Monde puzzle [#1051]

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

A combinatoric Le Monde mathematical puzzle of limited size:
When the only allowed move is to switch two balls from adjacent boxes, what is the minimal number of moves to return all balls in the above picture to their respective boxes? Same question with six boxes and 12 balls.

The question is rather interesting to code as I decided to use recursion (as usual!) but wanted to gain time by storing the number of steps needed by any configuration to reach its ordered recombination. Meaning I had to update an external vector within the recursive function for each new configuration I met. With help from Julien Stoehr, who presented me with the following code, a simplification of a common R function

v.assign <- function (i,value,...) {
  temp <- get(i, pos = 1)
  temp[...] <- value
  assign(i, temp, pos = 1)}

which assigns one or several entries to the external vector i. I thus used this trick in the following R code, where cosz is a vector of size 5¹⁰, much larger than the less than 10! values I need but easier to code. While n≤5.

n=5;tn=2*n
baz=n^(0:(tn-1))
cosz=rep(-1,n^tn)
swee <- function(balz){
  indz <- sum((balz-1)*baz)
  if (cosz[indz]==-1){ 
  if (min(diff(balz))==0){ #ordered
     v.assign("cosz",indz,value=1)}else{
       val <- n^tn
       for (i in 2:n)
       for (j in (2*i-1):(2*i))
       for (k in (2*i-3):(2*i-2)){
         calz <- balz
         calz[k] <- balz[j];calz[j] 0) 
           val <- min(val,1+swee(calz))}
     v.assign("cosz",indz,value=val)
  }}
 return(cosz[indz])}

which returns 2 for n=2, 6 for n=3, 11 for n=4, 15 for n=5. In the case n=6, I need a much better coding of the permutations of interest. Which is akin to ranking all words within a dictionary with letters (1,1,…,6,6). After some thinking (!) and searching, I came up with a new version, defining

parclass=rep(2,n)
rankum=function(confg){
    n=length(confg);permdex=1
    for (i in 1:(n-1)){
      x=confg[i]
      if (x>1){
        for (j in 1:(x-1)){
            if(parclass[j]>0){
                parclass[j]=parclass[j]-1
                permdex=permdex+ritpermz(n-i,parclass)
                parclass[j]=parclass[j]+1}}}
        parclass[x]=parclass[x]-1}
    return(permdex)}

ritpermz=function(n,parclass){
    return(factorial(n)/prod(factorial(parclass)))}

for finding the index of a given permutation, between 1 and (2n)!/2!..2!, and then calling the initial swee(p) with this modified allocation. The R code was still running when I posted this entry… and six days later, it returned the answer of 23.

Le Monde puzzle [#1043]

Posted in Books, Kids with tags , , , , , , on March 5, 2018 by xi'an

An arithmetic Le Monde mathematical puzzle :

A number is “noble” if all its digits are different and if it is equal to the average of all numbers created by permuting its digits. What are the noble numbers?

There is no need for simulation when plain enumeration works. After failing to install the R packge permutations, I installed the R package permute, which works, although (a) the function allPerm does not apply directly to a vector of characters or numbers but only to its size:

> allPerms(c("a","r","h"))
     [,1] [,2] [,3]
[1,]    1    3    2
[2,]    2    1    3
[3,]    2    3    1
[4,]    3    1    2
[5,]    3    2    1

and (b) as seen above the function does not contain “all” permutations since it misses the identity permutation.  Which ends up being fine for solving this puzzle. Using a bit of digit-character manipulation

findzol=function(N=2){
  for (u in 1:(10^N-1)){
    digz=strsplit(as.character(u),"")[[1]]
    if (length(digz)<N) 
      digz=c(rep("0",N-length(digz)),digz)
    if (length(unique(digz))==N){
      permz=apply(matrix(digz[allPerms(1:N)],
             ncol=N),2,as.numeric)
      if (mean(permz%*%10^{(N-1):0})==u) print(u)}}}

I found solutions for N=3

> findzol(3)
[1] 370
[1] 407
[1] 481
[1] 518
[1] 592
[1] 629

and none for N=4,5,6. Le Monde gives solutions for N=9, which is not achievable by my code!