G Fort

Member, IEEE Barbara Pascal
email: barbara.pascal@cnrs.fr

Fellow, IEEE Patrice Abry

Member, IEEE. Nelly Pustelnik

Covid19 Reproduction Number: Credibility Intervals by Blockwise Proximal Monte Carlo Samplers

Keywords: Markov Chain Monte Carlo sampling, nonsmooth convex optimization, Bayesian inverse problems, credibility intervals, Covid19, reproduction number

Monitoring the Covid19 pandemic constitutes a critical societal stake that received considerable research efforts. The intensity of the pandemic on a given territory is efficiently measured by the reproduction number, quantifying the rate of growth of daily new infections. Recently, estimates for the time evolution of the reproduction number were produced using an inverse problem formulation with a nonsmooth functional minimization. While it was designed to be robust to the limited quality of the Covid19 data (outliers, missing counts), the procedure lacks the ability to output credibility interval based estimates. This remains a severe limitation for practical use in actual pandemic monitoring by epidemiologists that the present work aims to overcome by use of Monte Carlo sampling. After interpretation of the nonsmooth functional into a Bayesian framework, several sampling schemes are tailored to adjust the nonsmooth nature of the resulting posterior distribution. The originality of the devised algorithms stems from combining a Langevin Monte Carlo sampling scheme with Proximal operators. Performance of the new algorithms in producing relevant credibility intervals for the reproduction number estimates and denoised counts are compared. Assessment is conducted on real daily new infection counts made available by the Johns Hopkins University. The interest of the devised monitoring tools are illustrated on Covid19 data from several different countries.

I. INTRODUCTION

Context. The Covid19 pandemic is causing unprecedented health, social, and economic crises. This triggered massive research efforts to design efficient procedures aiming to assess the intensity of the pandemic, a prerequisite to develop efficient sanitary policies [START_REF] Flahault | COVID-19 cacophony: is there any orchestra conductor?[END_REF]. Several indices are commonly used to measure the strength of a pandemic, such as, e.g., the reproduction number of interest here. However, often, the value of the index alone is not sufficient and credibility intervals of these indices constitute valuable information for the decision makers, notably in periods of rapid pandemic evolution or of changes in trends, an issue not always addressed in pandemic monitoring and at the heart of the present work.

Related works. Pandemic monitoring can be conducted with numerous tools from different scientific fields, (cf. [START_REF] Arino | Describing, Modelling and Forecasting the Spatial and Temporal Spread of COVID-19: A Short Review[END_REF] for a review), amongst which compartmental models, such as the founding Susceptible-Infectious-Recovered scheme. Within pandemic period, when data are scare and of limited quality, the reproduction number, R t , is often used by epidemiologists as an efficient practical proxy for the pandemic intensity: it measures the number of second infections caused by one primary infection (cf. e.g., [START_REF] Diekmann | On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations[END_REF][START_REF] Wallinga | Different Epidemic Curves for Severe Acute Respiratory Syndrome Reveal Similar Impacts of Control Measures[END_REF][START_REF] Van Den Driessche | Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF][START_REF] Obadia | The R0 package: A toolbox to estimate reproduction numbers for epidemic outbreaks[END_REF][START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF]). It thus plays a key role in the pandemic evolution assessment: the number of new infections today, Z t , depends on R t and on a weighted average of the new infection counts on previous days {. . . , Z t-3 , Z t-2 , Z t-1 }; the weight function being the so-called serial interval function Φ, that quantifies the distribution of the random delays between the onsets of symptoms in a primary and secondary cases [START_REF] Obadia | The R0 package: A toolbox to estimate reproduction numbers for epidemic outbreaks[END_REF][START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF][START_REF] Thompson | Improved inference of time-varying reproduction numbers during infectious disease outbreaks[END_REF][START_REF] Liu | Measurability of the epidemic reproduction number in data-driven contact networks[END_REF]. It has recently been shown that, within pandemic, reliable estimates for the temporal estimation of R t can be obtained from an inverse problem formulation resulting in a nonsmooth convex optimization problem [START_REF] Abry | Spatial and temporal regularization to estimate COVID-19 reproduction number R(t): Promoting piecewise smoothness via convex optimization[END_REF][START_REF] Pascal | Nonsmooth convex optimization to estimate the covid-19 reproduction number space-time evolution with robustness against low quality data[END_REF]. The functional to minimize is built from combining the pandemic model in [START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF], with time regularity constraints. While the procedure was engineered to produce realistic estimations of the temporal evolution of the reproduction number that are robust to the limited quality of the Covid19 pandemic data (severely corrupted with outliers, missing or negative counts and pseudo-seasonalities), it does not however provide credibility intervals, a critical issue towards its practical and actual use by epidemiologists, the main issue that we aim to address in the present work. Goals, contributions and outline. The overall goal of the present work is to devise Monte Carlo sampling strategies to perform the estimation by means of credibility intervals of the pandemic reproduction number and of denoised infection counts. To that end, Section II details the proposed statistical model used to embed into a stochastic framework the epidemiological model in [START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF] and its robust extension to data corruption [START_REF] Pascal | Nonsmooth convex optimization to estimate the covid-19 reproduction number space-time evolution with robustness against low quality data[END_REF]. Its originality stems from using non-differentiable priors to ensure robustness to data corruption. The uniqueness of the maximum a posteriori is thoroughly studied. Further, Section III devises original sampling schemes tailored to handle the non-differentiability of the target distribution. We propose two blockwise Proximal-Gradient based extensions of the Langevin Metropolis algorithms: PGdec and PGdual. We establish their ergodicity, and carry out a comparative study. Using real Covid19 data, made available at the Johns Hopkins University repository and described in Section IV, the performance of up to twelve variations of the sampling strategies are assessed and compared, using well-thought indices quantifying their efficiency (cf. Section V). Finally, in Section VI, the relevance of the proposed blockwise Proximal-Gradient samplers is illustrated for several different countries representative of the evolution of the pandemic across the world, for a 5week recent period. Daily updates of these credibility interval estimates as well as MATLAB routines for their calculations are available at https://github.com/gfort-lab/OpSiMorE. Notations. Vectors are column-vectors, and are denoted in bold font, e.g., R ∈ R T . For p ≤ q, the vector x p:q concatenates the scalars x i or vectors x i for i = p, . . . , q. The matrices are denoted in bold font. For a matrix A, A ⊤ (resp. det(A) and A -1 ) denotes the transpose of A (resp. the determinant and the inverse of A). We set A -⊤ := (A ⊤ ) -1 = (A -1 ) ⊤ . I p is the p × p identity matrix, and 0 p×q is the p × q null matrix (0 p×1 will be denoted by 0 p ). For a vector x ∈ R p , x 1 is the L 1 -norm and x is the L 2 -norm. Finally, N r (µ, C) denotes the R r -valued Gaussian distribution with expectation µ and covariance matrix C. For some γ > 0, the proximity operator of a proper, convex, lower semi-continuous function f from R d to ] -∞, +∞] is defined as

(∀x ∈ R d ) prox γf (x) := arg min y∈R d γf (y) + 1 2 y -x 2 .

II. COVID19 PANDEMIC BAYESIAN MODEL

A. Pandemic model

The present work makes use of a pandemic model devised by epidemiologists in [START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF] that focuses on a main pandemic index: the reproduction numbers, to be estimated from daily new infection counts. Elaborating on [START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF], it was further proposed in [START_REF] Pascal | Nonsmooth convex optimization to estimate the covid-19 reproduction number space-time evolution with robustness against low quality data[END_REF] to account for the limited quality of the intra-pandemic Covid19 data -highly corrupted by irrelevant, missing and mis-reported counts or by pseudo-seasonal effects -by means of additional outliers O, also unknown and to be estimated. The goal of the present work is thus to estimate, from a vector of T observed daily new infection counts

Z := (Z 1 , . . . , Z T ) ⊤ ∈ N T , the vector of unknowns (R, O) := (R 1 , . . . , R T ) ⊤ , (O 1 , . . . , O T ) ⊤ ∈ (R + ) T × R T
gathering the reproduction numbers R t and the outliers O t of T consecutive days.

B. A statistical model

Estimation entails the recourse to a formulation of the pandemic model, where the unknown parameters

θ := (R, O),
are a realization of a random vector; the statistical model defines its probability distribution. We consider a density with respect to (w.r.t.) the Lebesgue measure on (R

+ ) T × R T of the form θ → π Z (θ) ∝ exp (-f Z (θ) -g(θ)) ½ D Z (θ), (1)
where

½ A denotes the {0, 1}-valued indicator function of the set A. We choose g, f Z , D Z such that θ → -ln π Z (θ) on D Z
is the criterion proposed in [START_REF] Pascal | Nonsmooth convex optimization to estimate the covid-19 reproduction number space-time evolution with robustness against low quality data[END_REF] for the estimation of θ by minimization of a contrast. The function g is given by

g(θ) := λ R DR 1 + λ O O 1 ; (2)
D is proportional to the (T -2)×T discrete-time second order derivative matrix:

D := 1 √ 6     1 -2 1 0 0 . . . 0 0 1 -2 1 0 . . . 0 . . . . . . 0 . . . 1 -2 1     . (3)
The function f Z is given by

f Z (θ) := T t=1 {-Z t ln I t (θ, Z) + I t (θ, Z)} , (4)
where

I t (θ, Z) is a positive intensity. Up to an additive constant, f Z (θ) = t d KL (Z t |I t (θ, Z))
, where d KL denotes the Kullback-Leibler divergence related to the log-likelihood of a Poisson process 1 . We choose

I t (θ, Z) := R t τ φ u=1 Φ u Z t-u + O t . (5)
Φ := (Φ u ) 1≤u≤τ φ is the serial interval function, describing the average infectiousness profile after infection [START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF][START_REF] Thompson | Improved inference of time-varying reproduction numbers during infectious disease outbreaks[END_REF][START_REF] Liu | Measurability of the epidemic reproduction number in data-driven contact networks[END_REF]. Φ is assumed known and, following [START_REF] Riccardo | Epidemiological characteristics of COVID-19 cases in Italy and estimates of the reproductive numbers one month into the epidemic[END_REF][START_REF] Guzzetta | The impact of a nation-wide lockdown on COVID-19 transmissibility in Italy[END_REF], is classically modeled as a Gamma distribution truncated over τ φ = 26 days with mean and standard deviation of 6.6 and 3.5 days. The reproduction number R t at time #t acts on the rate at which a person infected at time t-u, generates new infections at time t: this rate is equal to R t Φ u . The support of the density π Z (•) is the measurable subset D Z of (R + ) T × R T ensuring that the intensities I t (θ, Z) are positive, or non-negative when Z t = 0 (with the convention that 0 ln 0 = 0 in the expression of f Z ). This yields

D Z := {θ : I t (θ, Z) > 0 for t s.t. Z t > 0} ∪ {θ : I t (θ, Z) ≥ 0 for t s.t. Z t = 0}. (6)

C. A Bayesian model

Let us describe a Bayesian framework such that π Z is a posterior distribution. The model depends on initial values Z 1-τ φ :0 , R 1 , R 2 , which are omitted in the notations. Prior distribution. Conditionally to (R 1 , R 2 ), the prior distribution of (R 3:T , O) has the density

T t=3 p 1 (R t |R t-1 , R t-2 ) T t=1 p 2 (O t )
1 Its definition is, for some z ∈ N,

(∀I ∈ R) dKL(z|I) :=      z ln z I + I -z if z > 0, I > 0, I if z = 0, I ≥ 0, +∞ otherwise.
w.r.t. the Lebesgue measure on R T -2 × R T ; we set

p 1 (r t |r t-1 , r t-2 ) := λ R 2 √ 6 exp - λ R √ 6 |r t -2r t-1 + r t-2 | p 2 (o t ) := λ O 2 exp(-λ O |o t |).
It is easily seen that this prior density is equal to exp(-g(θ)) up to a multiplicative constant (depending on λ R and λ O ). Such a prior states that R 3:T and O are mutually independent conditionally to (R 1 , R 2 ). Further, the outliers O are independent and distributed as a Laplace distribution (with parameter λ O > 0) as commonly encountered in the literature (see, e.g., [START_REF] Moulin | Analysis of multiresolution image denoising schemes using generalised Gaussian and complexity priors[END_REF][START_REF] Figueiredo | Adaptive sparseness for supervised learning[END_REF][START_REF] Park | The Bayesian Lasso[END_REF]). Finally R 3:T is distributed as an AR(2) process with Laplace distribution (with parameter λ R > 0). λ R and λ O are (fixed) positive regularization hyperparameters. g favors some large values among many small ones of the components of both vectors DR and O; it models smooth piecewise linear time evolutions for R, or equivalently a sparse set of components where the discrete second order derivative in time of R is non zero.

Likelihood. Conditionally to (R, O), the observations are not independent and the joint distribution of Z is

T t=1 p 3 (Z t |Z t-τ φ :t-1 , R t , O t ).
When R t ≥ 0 and I t (θ, Z) > 0 (see [START_REF] Van Den Driessche | Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF] for the definition of I t ), following [START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF][START_REF] Pascal | Nonsmooth convex optimization to estimate the covid-19 reproduction number space-time evolution with robustness against low quality data[END_REF], the pandemic diffusion is modeled as a Poisson distribution with time-varying intensity I t (θ, Z):

p 3 (Z t |Z 1:t-1 , R t , O t ) := I t (θ, Z) Zt Z t ! exp(-I t (θ, Z)). ( 7
)
Such a model claims that the mean value of the counts Z t at time t is

τ φ u=1 R t Φ u Z t-u +O t .
When R t ≥ 0 and I t (θ, Z) = 0, p 3 is the Dirac mass at zero: it is equal to [START_REF] Cori | A new framework and software to estimate time-varying reproduction numbers during epidemics[END_REF] when Z t = 0 and equal to 0 otherwise. Finally, when R t < 0 or I t (θ, Z) < 0, p 3 is a distribution on the negative integers; this implies that when

Z ∈ N T , p 3 (Z t |Z t-τ φ :t-1 , R t , O t ) = 0.
This description implies that the likelihood of a vector of observations Z taking values in

N T is exp(-f Z (θ))½ D Z (θ). Posterior distribution. This model implies that π Z (θ) is the a posteriori distribution of (R 3:T , O) given Z when Z ∈ N T .

D. Bayesian estimators

In the Bayesian approach to Decision Theory, the maximum, the median, and the expectation of the a posteriori distribution, are Bayes estimators θ associated to a loss function ℓ

θ(Z) := Argmin τ ∈D Z D Z ℓ(τ , θ)π Z (θ)dθ;
ℓ is, respectively, the 0 -1 loss, the L 1 -norm and the squared L 2 -norm (see e.g. [START_REF] Robert | The Bayesian choice: a decision-theoretic motivation[END_REF]Sections 2.3. and 2.5.]). Computing the Maximum a Posteriori (MAP) fits the minimization problem proposed in [START_REF] Pascal | Nonsmooth convex optimization to estimate the covid-19 reproduction number space-time evolution with robustness against low quality data[END_REF] for the reconstruction of θ:

Argmin R,O T t=1 d KL (Z t |I t (θ, Z)) + λ R DR 1 + λ O O 1 . (8)
The optimization problem is a nonsmooth convex minimization problem both encapsulating the transmission process, favoring piecewise linear behavior of the reproduction number along time and sparsity of the outliers. The minimization is performed with the Chambolle-Pock primal-dual algorithm allowing to handle both the non-differentiability and linear operators [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. Properties of the MAP are established in Proposition 1.

Proposition 1. If there are at least two positive averaged counts

τ φ u=1 Φ u Z t⋆-u , τ φ u=1 Φ u Z t⋆⋆-u , and one positive count Z τ , a MAP exists. If θ ⋆ = (R ⋆ , O ⋆ ) and θ ⋆⋆ = (R ⋆⋆ , O ⋆⋆ ) are two MAP, then I t (θ ⋆ , Z) = I t (θ ⋆⋆ , Z), O ⋆
t O ⋆⋆ t ≥ 0 and (DR ⋆ ) t (DR ⋆⋆ ) t ≥ 0 for any t ∈ {1, . . . , T }. Proof. The first statement is adapted from [START_REF] Pascal | Nonsmooth convex optimization to estimate the covid-19 reproduction number space-time evolution with robustness against low quality data[END_REF]; The second statement is established in [START_REF] Pascal | Nonsmooth convex optimization to estimate the covid-19 reproduction number space-time evolution with robustness against low quality data[END_REF]. The sign conditions result from a first order expansion of the L 1 -norm. For a detailed proof, see section VIII in the Supplementary material.

Proposition 1 implies that the MAP is either unique, or that there are uncountably many MAP. In addition, it shows that f Z and g are constant over the set of the minimizers. Thus, following the same lines as in [START_REF] Ali | The Generalized Lasso Problem and Uniqueness[END_REF], a sufficient condition for the uniqueness of the MAP is derived (see section VIII in the Supplementary material).

The expression of the distribution π Z in ( 1) is so complex that it is known only up to a normalizing constant. Consequently, the computation of most statistics of π Z (•) relies on Monte Carlo samplers, in order to produce samples {θ n , n ≥ 0} in D Z approximating π Z (see e.g. [21, section 2.3]): for example, the estimation of the median and more generally the quantiles of π Z can rely on the order statistics of the samples, and the mean a posteriori can be approximated by the Monte Carlo sum N -1 N n=1 θ n .

III. BLOCKWISE PROXIMAL-GRADIENT MONTE CARLO

SAMPLERS

The aim is now to devise Monte Carlo sampling strategies for the distribution defined in [START_REF] Flahault | COVID-19 cacophony: is there any orchestra conductor?[END_REF]. However, this section will address a broader class of densities, defined on R d with respect to the Lebesgue measure, and expressed as π(θ) ∝ exp(-F (θ))½ D (θ) where F := f + g and f , g, and D satisfy the smoothness and blockwise structure assumptions A1 and A2 defined below.

A. Smoothness and Blockwise structure

A 1. f and g are finite on D ⊆ R d and f is continuously differentiable on the interior of D.

Additionally, g has a blockwise structure that we aim to use in the design of the proposed samplers. This blockwise structure stems both from the decomposition of θ into J blocks (θ 1 , . . . , θ J ) ∈ R d1 × . . . × R dJ and from the sum of several functions of θ j possibly combined with a linear operator.

A2. For j ∈ {1, . . . , J}, i ∈ {1, . . . , I j }, there exist matrices A i,j ∈ R ci,j ×dj , and proper, convex, lower semi-continuous functions g i,j : R ci,j → ] -∞, +∞] such that

J j=1 d j = d and ∀θ := (θ ⊤ 1 , . . . , θ ⊤ J ) ⊤ , g(θ) := J j=1 Ij i=1
g i,j (A i,j θ j ) .

In addition, the proximity operator of g i,j has a closed form expression.

In Bayesian inverse problems, f may stand for the data fidelity term and the nonsmooth part g stands for many penalty terms acting on blocks of the parameter θ. Different splittings of the prior defined by (2) fits Assumption A2.

Example 2. The prior g given by (2) satisfies A2:

λ R DR 1 + λ O O 1 = g 1,1 (A 1,1 R) + g 1,2 (A 1,2 O),
where

θ 1 := R, θ 2 := O, A 1,1 := D, A 1,2 = I T , g 1,1 := λ R • 1 , and g 1,2 := λ O • 1 .
Example 3. The prior g given by (2) satisfies A2:

λ R DR 1 + λ O O 1 = 3 i=1 g i,1 (A i,1 R) + g 1,2 (A 1,2 O),
where θ 1 := R and θ 2 := O. A i,1 collects the rows i, i + 3, i + 6, . . . , of the matrix D, A 1,2

:= I T , g i,1 := λ R • 1 and g 1,2 := λ O • 1 .
The second example follows block splitting strategies described in [START_REF] Pustelnik | Parallel proXimal algorithm for image restoration using hybrid regularization[END_REF][START_REF] Pascal | Block-coordinate proximal algorithms for scale-free texture segmentation[END_REF].

B. Proximal algorithms and Metropolis-Hastings algorithms

The design of an optimization strategy to minimize F on D and the design of a sampler to approximate the target distribution π both rely on the activation of an operator µ : R d → R d . To be more specific, when minimizing F , we aim to design a sequence of the form:

θ n+1 = µ(θ n ) ( 9
)
where µ is an operator built from F and D in such a way that the sequence (θ n ) n∈N converges to a minimizer of F (cf. [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Chambolle | An introduction to continuous optimization for imaging[END_REF] for an exhaustive list of algorithmic schemes). When building a Metropolis-Hastings algorithm (say with Gaussian proposal), a new point is proposed as:

θ n+1/2 = µ(θ n ) + ξ n+1 where ξ n+1 ∼ N d (0 d , C); (10
) C ∈ R d×d is a positive definite matrix. The Langevin dynamics is recovered in the specific case where F is smooth with µ(θ) := θ -γ∇F (θ) being a gradient ascent over ln π with step size γ > 0 and C := 2γI d . Scaled Langevin samplers are also popular: given a d × d matrix Γ, set

µ(θ) := θ -γΓΓ ⊤ ∇F (θ), C := 2γ ΓΓ ⊤ ; (11)
they are inherited from the so-called tempered Langevin diffusions [START_REF] Kent | Time-Reversible Diffusions[END_REF] (see also [START_REF] Roberts | Langevin Diffusions and Metropolis-Hastings Algorithms[END_REF] for a pioneering work on its use in the Markov Chain Monte Carlo literature). Either this proposed point is the new point θ n+1 = θ n+1/2 (thus yielding the Langevin Monte Carlo algorithm [START_REF] Parisi | Correlation functions and computer simulations[END_REF]; see also [START_REF] Durmus | Nonasymptotic convergence analysis for the unadjusted Langevin algorithm[END_REF][START_REF] Dwivedi | Logconcave sampling: Metropolis-Hastings algorithms are fast![END_REF]) or there is an acceptance-rejection Metropolis mechanism (thus yielding the Metropolis Adjusted Langevin Algorithm (MALA) [START_REF] Roberts | Exponential convergence of Langevin distributions and their discrete approximations[END_REF]). The general Metropolis-Hastings procedure with Gaussian proposal, is summarized in Algorithm 1 where we denote by q(θ, τ ) the density of the distribution

N d (µ(θ), C) evalu- ated at τ ∈ R d : q(θ, τ ) := exp -0.5 (τ -µ(θ)) ⊤ C -1 (τ -µ(θ)) √ 2π d det(C)
.

The constraint θ ∈ D is managed by the acceptance-rejection step (since π(θ n+1/2 ) = 0 when θ n+1/2 / ∈ D) but not necessarily in the proposal mechanism. The MALA algorithms Algorithm 1: Metropolis-Hastings with Gaussian proposal Data: a positive definite matrix C; γ > 0; a positive integer

N max ; θ 0 ∈ D Result: A D-valued sequence {θ n , n ∈ {0, . . . , N max }} 1 for n = 0, . . . , N max -1 do 2 Sample ξ n+1 ∼ N d (0 d , C); 3 Set θ n+ 1 2 = µ(θ n ) + ξ n+1 ; 4 Set θ n+1 = θ n+ 1 2 with probability 1 ∧ π(θ n+ 1 2 ) π(θ n ) q(θ n+ 1 2 , θ n ) q(θ n , θ n+ 1 2 ) (12)
and θ n+1 = θ n otherwise.

drift the proposed moves towards areas of high probability for the distribution π, using first order information on π. Building on this idea, many strategies were proposed in the literature in the setting defined by A1: µ can either be a gradient step when F is smooth, or a proximal step (i.e., µ = prox γF also referred as an implicit subgradient descent step), or a Moreau-Yosida envelope gradient step (i.e. µ = I d -γ∇( γ F ) where the Moreau envelope of a function F with parameter γ > 0 is defined as γ F := inf y γF (y) + 1 2 • -y 2 ). In [START_REF] Durmus | Analysis of Langevin Monte Carlo via Convex Optimization[END_REF], explicit subgradient steps possibly combined with a proximal step are used. In [START_REF] Chatterji | Langevin Monte Carlo without smoothness[END_REF], µ relies on a Gaussian smoothing of convex functions with Hölder-continuous sub-gradients; this method applies under regularity conditions and convexity assumptions on g which are not implied by A1-A2. In [START_REF] Durmus | Efficient Bayesian Computation by Proximal Markov Chain Monte Carlo: When Langevin Meets Moreau[END_REF], the authors add a Moreau-Yosida envelope term and a gradient term. In [START_REF] Luu | Sampling from Non-smooth Distributions Through Langevin Diffusion[END_REF], µ composes a Moreau-Yosida envelope of g(A•) and a gradient step. Let us cite [START_REF] Atchadé | A Moreau-Yosida approximation scheme for a class of high-dimensional posterior distributions[END_REF][START_REF] Schreck | A Shrinkage-Thresholding Metropolis Adjusted Langevin Algorithm for Bayesian Variable Selection[END_REF][START_REF] Salim | Primal Dual Interpretation of the Proximal Stochastic Gradient Langevin Algorithm[END_REF] who also use proximal operators in order to define trans-dimensional Monte Carlo samplers -an objective which is out of the scope defined by A1. See also [START_REF] Pereyra | A survey of stochastic simulation and optimization methods in signal processing[END_REF] and [41, section 3] for a survey on MCMC samplers using optimization techniques, and [START_REF] Liu | The Multiple-Try Method and Local Optimization in Metropolis Sampling[END_REF]Section 4] for optimization techniques combined with Multiple Try Metropolis strategies.

However, in the optimization context ( 9) with a nonsmooth objective function F , explicit sub-gradient method is suboptimal with rate O(1/ √ n) compared to proximal-based strategy such as forward-backward (relying on an implicit subgradient method) whose rate is O(1/n) [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF]. Proximal-based algorithms are thus especially adapted to handle non-smooth functions but also smooth functions that do not have a Lipschitz gradient. The only limitation of using proximal-based algorithms relies on the possible difficulty to derive a closedform expression of the proximity operator. However, for most common functions encountered in signal processing, closedforms exist and have been summarized on the website Prox-Repository http://proximity-operator.net/scalarfunctions.htmlsee also Lemma 4 below, for standard functions composed with linear operators. When F = f + g where f and g satisfy Assumption A1, deriving a proximal activation µ is often a tedious task as no closed form expression for prox f +g exists in a general framework [START_REF] Yu | On decomposing the proximal map[END_REF][START_REF] Pustelnik | Proximity operator of a sum of functions; application to depth map estimation[END_REF]. The standard solution consists in a proximal-gradient activation: µ(θ) := prox γg (θ -γ∇f (θ)). When dealing with a blockwise structure for g as in A2, the choice of µ has to manage both the additive structure of g and the combination of the g i,j 's with a linear operator. Unfortunately, the proximity operator has a closed form expression in very limited cases recalled below in Lemma 4.

Lemma 4. Let A ∈ R c×d . 1) Let g := 1 2 • -z 2 with z ∈ R c . For every γ > 0, prox γg(A•) = (γA ⊤ A + I d ) -1 (• + γA ⊤ z).
2) Let g be a proper lower semi-continuous convex function. Let AA ⊤ be invertible. For every γ > 0,

prox γg(A•) = I d -A ⊤ (AA ⊤ ) -1 (I d -prox (AA ⊤ ) -1 γg )A
where, for every

x ∈ R d , prox (AA ⊤ ) -1 γg (x):=argmin y g(y)+ 1 2γ (x-y) ⊤ (AA ⊤ ) -1 (x-y).
3) Let g(A•) := c ℓ=1 g ℓ (a ℓ •) where g ℓ is convex, lower semi-continuous, and proper from R to ] -∞, +∞], and a ℓ ∈ R d denotes the row #ℓ of A. Suppose that AA ⊤ = Λ, where Λ := diag(χ 1 , . . . , χ c ) and χ ℓ > 0. Then, for every γ > 0:

for all η ∈ R d prox γg(A•) (η) = η -A ⊤ Λ -1 Aη -prox γΛg (Aη) where for all ζ = ζ 1:c , we set Λg(ζ) := c ℓ=1 χ ℓ g ℓ (ζ ℓ ).
Proof. 1) See [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. 2) This result is a direct consequence of [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Proposition 23.25 (ii)] and [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Proposition 23.345 [START_REF] Pustelnik | Parallel proXimal algorithm for image restoration using hybrid regularization[END_REF] and a direct consequence of 2) for specific choices of A and g.

(ii)- (iii)]. 3) Result extracted from

C. Novel Blockwise Proximal Monte Carlo Samplers: PGdec and PGdual

None of the algorithms recalled in the previous section directly applies to the context of Assumptions A1 to A2. We propose two novel Metropolis-Hastings algorithms: the PGdec sampler and the PGdual sampler, which use the proximal operator to handle the nonsmooth part g of ln π, its additive structure and its combination with linear operators. Such strategies can be used as building blocks for more general MCMC samplers such as Multiple Try based samplers (see e.g. [START_REF] Martino | A review of multiple try MCMC algorithms for signal processing[END_REF] and [START_REF] Bédard | Scaling analysis of multiple-try MCMC methods[END_REF]) and more generally, adaptive Monte Carlo samplers (see e.g. [START_REF] Andrieu | A tutorial on adaptive MCMC[END_REF]); An extension to Gibbs sampler is proposed and discussed in Section III-F.

• The PGdec sampler. Additionally to Assumptions A1 and A2, we further assume: A3. Each function g i,j (A i,j •) possesses a proximal operator having a closed form expression.

A3 assumes that each component g i,j (A i,j •) has a tractable proximal operator which does not imply that the sum g = J j=1 Ij i=1 g i,j admits a tractable proximal operator. The Proximal-Gradient Decomposition sampler (PGdec) is described by algorithm 2. It is a Metropolis-Hastings sampler with Gaussian proposal: conditionally to the current point θ n , PGdec proposes a move to θ n+1/2 = (θ

n+1/2 1 , • • • , θ n+1/2 J
) sampling independently the J blocks from Gaussian distributions (see θ n+1/2 j in line 5 of algorithm 2); Then an acceptance-rejection step is applied (see [START_REF] Figueiredo | Adaptive sparseness for supervised learning[END_REF]). The originality of our method is the definition of µ: for every j ∈ {1, . . . , J} and i ∈ {1, . . . , I j },

µ PGdec i,j (θ) := prox γj gi,j (Ai,j •) (θ j -γ j ∇ j f (θ)) , (13)
where γ j is a positive step size and ∇ j denotes the differential operator w.r.t. the block #j of θ. The proposed drift takes benefit of the blockwise separable expression of g and computes at each iteration the proximal operator associated to part of the sum in order to perform the proximal activation with a closed form expression. Conditionally to θ n , for each block #j, one of the component #i j is selected at random in {1, . . . , I j } (see line 3); then, θ n+1/2 j is sampled from a R djvalued Gaussian distribution with expectation µ PGdec ij ,j (θ n ) and covariance matrix C ij ,j . We denote by q i,j (θ, τ ) the density of the distribution N dj (µ i,j (θ), C i,j ) evaluated at τ ∈ R dj :

q i,j (θ, τ ) := exp -0.5 (τ -µ i,j (θ)) ⊤ C -1 i,j (τ -µ i,j (θ)) √ 2π dj det(C i,j )
.

Remark 5. Let j ∈ {1, . . . , J} and i ∈ {1, . . . , I j }. Given

θ n = (θ n 1:J ) ∈ R d , µ PGdec i,j
(θ n ) successively computes a gradient step w.r.t. the smooth function f and the variable θ j , and a proximal step with respect to the function g i,j (A i,j •); the step size is γ j for both steps. Hence, µ PGdec i,j

(θ n ) is a Proximal- Gradient (PG) step w.r.t. the function θ j → f (θ n 1:j-1 , θ j , θ n j+1:J ) + g i,j (A i,j θ j ). ( 14
)
Example 6 (Example 2 to follow). The proximal operator of R → λ R DR 1 is not explicit, so that, when decomposing g as proposed in Example 2, PGdec can not be applied to approximate π Z .

Example 7 (Example 3 to follow). For all i ∈ {1, 2, 3}, A i,1 A ⊤ i,1 is the identity matrix so that, by Lemma 4, prox γ1 λ R Ai,1• 1 is explicit and given by (

I T -A ⊤ i,1 A i,1 ) + A ⊤ i,1 prox γ1 λ R • 1 (A i,1 •).
In addition, prox γ2 λ O • 1 has a closed form expression. Hence, when decomposing g as proposed in Example 3, PGdec can be applied to approximate π Z .

Algorithm 2: Blockwise Metropolis-Hastings samplers Data: d j × d j positive definite matrices C i,j ; γ j > 0; a positive integer N max ; θ 0 ∈ D Result: A D-valued sequence {θ n , n ∈ {0, . . . , N max }} 1 for n = 0, . . . , N max -1 do 2 for j = 1, . . . , J do 3 Sample i j ∈ {1, . . . , I j } with probability 1/I j ;

4 Sample ξ n+1 j ∼ N dj (0 dj , C ij ,j ); 5 Set θ n+ 1 2 j = µ ij ,j (θ n ) + ξ n+1 j ; 6 Set θ n+1 = θ n+ 1 2 with probability 1 ∧ π(θ n+ 1 2 ) π(θ n ) J j=1 q ij ,j (θ n+ 1 2 , θ n j ) q ij ,j (θ n , θ n+ 1 2 j ) (15)
and θ n+1 = θ n otherwise.

• The PGdual sampler. The Proximal-Gradient dual sampler (PGdual) is defined along the same lines as algorithm 2. It is designed for situations when for any i, j, the dimensions of the matrices A i,j satisfy c i,j ≤ d j and A i,j can be augmented in an invertible d j × d j matrix -denoted by Āi,j .

For every j ∈ {1, . . . , J} and i ∈ {1, . . . , I j }, let Āi,j be a d j × d j invertible matrix such that for any θ j ∈ R dj , ( Āi,j θ j ) dj -ci,j+1:dj = A i,j θ j ; For x = (x 1 , . . . , x dj ) ∈ R dj , define ḡi,j (x) := g i,j (x dj-ci,j +1:dj ). PGdual uses the drift functions

µ PGdual i,j (θ) := Ā-1 i,j prox γj ḡi,j Āi,j θ j -γ j Ā-⊤ i,j ∇ j f (θ) . ( 16
)
Remark 8. For every j ∈ {1, . . . , J}, select i j ∈ {1, . . . , I j }, and consider the partial objective function:

θ → f (θ 1 , . . . , θ J ) + j ḡij,j ( Āij,j θ j ). ( 17
)
Consider the one-to-one maps θ j := Āij,j θ j for any j, and the application

θ = θ 1:J → f ( Ā-1 i1,1 θ 1 , . . . , Ā-1 iJ ,J θ J ) + j ḡij,j ( θ j ). ( 18
)
The PG step w.r.t. θ j reads:

prox γj ḡi j ,j θ j -γ j Ā-⊤ ij ,j ∇ j f ( Ā-1 i1,1 θ 1 , . . . , Ā-1 iJ ,J θ J ) = prox γj ḡi j ,j Āij,j θ j -γ j Ā-⊤ ij ,j ∇ j f (θ) . ( 19
)
Therefore, applying a PG step w.r.t. the variable θ j in this "dual" space, and going back to the original space by applying Ā-1 ij ,j leads to the drift µ PGdual ij ,j in [START_REF] Park | The Bayesian Lasso[END_REF]. From [START_REF] Park | The Bayesian Lasso[END_REF] and Lemma 4 Item 2, we have

µ PGdual i,j (θ) := prox Ā⊤ i,j Āi,j γj ḡi,j ( Āi,j •) θ j -γ j ( Ā⊤ i,j Āi,j ) -1 ∇ j f (θ) (20)
thus showing that it is a Variable Metric Proximal-Gradient step in the metric induced by Ā⊤ i,j Āi,j , for the minimization of θ j → f (θ) + ḡi,j ( Āi,j θ j ). See [START_REF] Chen | Convergence rates in forward-backward splitting[END_REF][START_REF] Combettes | Variable metric forward-backward splitting with applications to monotone inclusions in duality[END_REF]; see also [51, section 2] in the specific case J = 1, I = 1.

Example 9 (Example 2 to follow). Denote by D I the T × T invertible matrix obtained by augmenting D with two vectors in R T as follows:

D I :=   1 0 0 . . . 0 -2/ √ 5 1/ √ 5 0 . . . 0 D   ; (21)
D is subdiagonal with normalized rows (see Eq. ( 3)) and we define D I so that it has the same properties. Then (D I R) 3:T = DR, and for

R ∈ R T , ḡ1,1 (R) = λ R R 3:T 1 . In addition, ḡ1,2 = g 1,2
. Hence, when decomposing g as proposed in Example 2, PGdual can be used to approximate π Z .

• Choice of the covariance matrix C. Based on Remark 8, Āi,j µ PGdual i,j

( θ) is a PG step in a dual space. Since such a step can be seen as an extension of a gradient step to a nonsmooth function, a natural idea is to mime the MALA proposal and add a Gaussian noise ξ with covariance matrix 2γ j I dj in the dual space. Therefore, in the original space, the noise ξ := Ā-1 i,j ξ has covariance matrix C i,j := 2γ j Ā-1 i,j

Ā-⊤ i,j . In ( 16), note the preconditioning matrix Ā-⊤ i,j before the gradient and the matrix Ā-1 i,j before the proximal operator: there is a parallel between such a choice of C i,j and the scaled MALA proposal mechanism [START_REF] Pascal | Nonsmooth convex optimization to estimate the covid-19 reproduction number space-time evolution with robustness against low quality data[END_REF]. Equivalently, note the parallel between the preconditioning of the gradient term in [START_REF] Ali | The Generalized Lasso Problem and Uniqueness[END_REF], the covariance matrix C i,j and ( 11).

An equivalent argumentation can be developed for µ PGdec .

D. Interpretations of PGdec and PGdual

In the case A i,j A ⊤ i,j = ν i,j I ci,j with ν i,j > 0, we have from Lemma 4 Item 3,

µ PGdec i,j (θ n ) = I dj -Π i,j G j (θ n ) + A ⊤ i,j A i,j A ⊤ i,j -1 prox νi,j γj gi,j (A i,j G j (θ n )) , (22)
where

G j (θ n ) := θ n j -γ j ∇ j f (θ n ) is a gradient step, and
Π i,j := A ⊤ i,j A i,j A ⊤ i,j -1 A i,j
is the orthogonal projection matrix on the range of A ⊤ i,j . [START_REF] Pustelnik | Parallel proXimal algorithm for image restoration using hybrid regularization[END_REF] shows that µ PGdec i,j (θ n ) is the sum of two orthogonal terms: the first one is the orthogonal projection of G j (θ n ) on the orthogonal space of the range of A ⊤ i,j , and the second term is in the range space of A ⊤ i,j . This second term may be seen as a proximal-contraction of Π i,j G j (θ n ).

Theorem 10 makes the PGdual idea explicit by proposing a matrix Āi,j augmenting A i,j , computing the associated drift µ PGdual i,j and comparing it to µ PGdec i,j . Theorem 10. Assume A1 and A2. Let j ∈ {1, . . . , J} and i ∈ {1, . . . , I j }. Assume that c i,j < d j and A i,j A ⊤ i,j is invertible. Let U i,j be a (d jc i,j ) × d j matrix such that U i,j A ⊤ i,j = 0 (dj -ci,j )×ci,j and U i,j U ⊤ i,j is invertible.

Then, the matrix Āi,j := [U i,j ; A i,j ] ∈ R dj ×dj is invertible. For any θ ∈ D, µ PGdual i,j

(θ) given by ( 16) is equal to

A ⊤ i,j A i,j A ⊤ i,j -1 prox γj gi,j A i,j θ j -γ j Ω i,j ∇ j f (θ) + (I dj -Π i,j ) (θ j -γ j Ω i,j ∇ j f (θ)) , (23)
where

Ω i,j := U ⊤ i,j (U i,j U ⊤ i,j ) -2 U i,j , Ω i,j := A ⊤ i,j (A i,j A ⊤ i,j ) -2 A i,j . Additionally, when U i,j U ⊤ i,j = I dj-ci,j , A i,j A ⊤ i,j = I ci,j , and under A3 then µ PGdual i,j = µ PGdec i,j
where µ PGdec i,j is given by [START_REF] Pustelnik | Parallel proXimal algorithm for image restoration using hybrid regularization[END_REF].

Proof. The main ingredients are the equalities

Ā-1 i,j = U ⊤ i,j (U i,j U ⊤ i,j ) -1 A ⊤ i,j (A i,j A ⊤ i,j ) -1 ;
and

U ⊤ i,j (U i,j U ⊤ i,j ) -1 U i,j + Π i,j = I dj .
The proof follows from standard matrix algebra. A detailed proof is given in section IX of the Supplementary material.

The result remains true when c i,j = d j ; in that case, Āi,j = A i,j , I dj -Π i,j = 0 dj×dj and Ω i,j = A -1 i,j A -⊤ i,j . Theorem 10 provides sufficient conditions for the drift function µ PGdec i,j and the drift function µ PGdual i,j to be equal. Observe that the conditions on U i,j are satisfied as soon as the rows of U i,j are orthonormal and orthogonal to the rows of A i,j .

Let us derive two strategies for the application of PGdual to sample the target density defined by [START_REF] Flahault | COVID-19 cacophony: is there any orchestra conductor?[END_REF].

D O := U 1,1 D ∈ R T ×T ( 24
)
where U 1,1 is obtained by making orthogonal the first two rows of D I (see [START_REF] Douc | Markov Chains, ser. Springer Series in Operations Research and Financial Engineering[END_REF]) and making them orthogonal to the rows of D. With this choice, D O has a condition number one order of magnitude smaller compared to D I , which is likely to impact the convergence of Markov chains (see Section VI). This strategy acts globally on ln π Z by proposing, at each iteration, a PG approach for the function

f Z + λ R D • 1 + λ O • 1 .
It will be numerically explored in section V.

Example 12 (Example 3 to follow).

A second strategy is to decompose g as in Example 3, and define the matrices Āi,j as described in Theorem 10. This strategy defines µ PGdual by considering part of ln π Z : at each iteration, having selected i 1 ∈ {1, 2, 3}, it proposes a PG approach for the function

f Z + λ R A i1,1 • 1 + λ O • 1 (see Remark 8).
E. Convergence analysis of PGdec and PGdual.

We prove that both PGdec and PGdual produce a sequence of points {θ n , n ≥ 0} which is an ergodic Markov chain having π as its unique invariant distribution. Remark 14. A slight adaption of the proof shows that when selecting the indices i j , the uniform distribution on {1, . . . , I j } can be replaced with a probability distribution {ρ i,j (θ n ), i = 1, . . . , I j } depending on the current value θ n , up to a modification of the acceptance-rejection ratio [START_REF] Figueiredo | Adaptive sparseness for supervised learning[END_REF] (see section IX of the Supplementary material).

F. Gibbs PGdec and Gibbs PGdual samplers.

The additive structure of g in A 2, naturally suggests the extension of PGdec and PGdual to the Gibbs sampler algorithm [START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF][START_REF] Gelfand | Sampling-based approaches to calculating marginal densities[END_REF]. Since whatever θ = θ 1:J ∈ R d , exact sampling from the conditional distributions on R dj τ → π j (τ |θ) ∝ exp(-f (θ 1:j-1 , τ , θ j+1:J ) -Ij i=1 g i,j (A i,j τ ))½ D (θ 1:j-1 , τ , θ j+1:J ) for j = 1, . . . , J, is not always explicit, we propose a Metropolis-within-Gibbs strategy; see [START_REF] Chan | Discussion: Markov Chains for Exploring Posterior Distributions[END_REF]. The pseudo-code of the so-called Gibbs Blockwise Proximal sampler is given by algorithm 3 in the case of a systematic scan order of the J components (our method easily extends to other scan orders; details are left to the reader): at each iteration #(n+1), and for each block #j, (i) sample at random i j ∈ {1, . . . , I j }, (ii) sample a candidate from N dj (µ ij ,j (ϑ), C ij ,j ) where ϑ := (θ n+1 1:j-1 , θ n j:J ) ∈ R d is the current value of the chain, and (iii) accept/reject this candidate via a Metropolis step targeting the distribution π j (•|ϑ).

The Gibbs PGdec and the Gibbs PGdual samplers correspond to algorithm 3 applied resp. with µ i,j = µ PGdec i,j and µ i,j = µ PGdual i,j .

IV. COVID-19 DATA

To illustrate, assess, and compare the relevance and performance of the Monte Carlo procedures for credibility interval estimation proposed here, use is made of real Covid-19 data made available by the Johns Hopkins University2 . The repository 3 collects, impressively since the outbreak of the pandemic,

π j (θ n j |ϑ) q ij ,j ((ϑ 1:j-1 , θ n+ 1 2 j
, ϑ j+1:J ), θ n j ) q ij ,j (ϑ, θ As mentioned in the Introduction section, because of the sanitary crisis context, data made available by most National Public Health Authorities in the world are of very limited quality as they are corrupted by outliers and missing or negative counts. Data quality also varies a lot across countries or even within a given country depending on the phases and periods of the pandemic.

The present work uses the new infection counts only, as the estimation of the space-time evolution of the pandemic reproduction number R is targeted.

A mild and non-informative preprocessing is applied to data by replacing negative counts by a null value.

V. MONTE CARLO SAMPLER ASSESSMENTS

This section aims to assess and compare the performance for several variations of the samplers PGdec and PGdual introduced in section III, using the real Covid19 data described in section IV. The target distribution is π Z , given by (1).

A. Assessment set-up

Performance assessments were conducted on data from several countries and for different periods of interests. For space reasons, they are reported only for the United Kingdom and for a recent time period (Dec. 6th, 2021 to Jan. 9th, 2022), with corresponding Z in Fig. 2[top right plot].

Following [START_REF] Pascal | Nonsmooth convex optimization to estimate the covid-19 reproduction number space-time evolution with robustness against low quality data[END_REF], we set λ O = 0.05 and λ R = 3.5 × √ 6 σ Z /4 with σ Z the standard deviation of Z. Samplers are run with N max = 1e7 iterations, including a burn-in phase of 3e6 iterations. Except for the plots in Fig. 1[first row], reported performance are computed by discarding the points produced during the burn-in phase. The initial values Z -τ φ +1 , . . . , Z -1 , Z 0 are set to the observed counts. For each algorithm, performances are computed from averages over 15 independent runs. Initial points consist of random perturbations around the non-informative point R 0 := (1, . . . , 1) ⊤ ∈ R T and O 0 := (0, . . . , 0) ⊤ ∈ R T . For all samplers, θ is seen as a two-block vectors (θ 1 = R, θ 2 = O); hence J = 2. For PGdec, g is decomposed as in Example 3; for PGdual, it is decomposed as in Example 2.

B. Samplers

PGdec and Gibbs PGdec. These two samplers are run by decomposing g as described in Example 3. This yields two different drift functions for the blocks R and O; see Example 7. For the R-part, the drift functions are defined, for i = 1, 2, 3, as:

µ PGdec i,1 (θ) := I T -A ⊤ i,1 A i,1 (R -γ 1 ∇ 1 f Z (θ)) + A ⊤ i,1 prox γ1λ R • 1 (A i,1 R -γ 1 A i,1 ∇ 1 f Z (θ)) , ( 25
)
obtained from [START_REF] Pustelnik | Parallel proXimal algorithm for image restoration using hybrid regularization[END_REF] with A i,1 A ⊤ i,1 the identity matrix. For the O-part,

µ PGdec 1,2 (θ) := prox γ2λ O • 1 (O -γ 2 ∇ 2 f Z (θ)) . (26)
∇ 1 (resp. ∇ 2 ) denotes the gradient operator w.r.t R (resp. O).

PGdual and Gibbs PGdual. These two samplers are run by decomposing g as described in Example 2. Consequently, for the R-part, the drift function is given by

µ PGdual 1,1 (θ) := D -1 prox γ1λ R (•)3:T 1 DR -γ 1 D -⊤ ∇ 1 f Z (θ) (27
µ PGdual 1,2 (θ) := prox γ2λ O • 1 (O -γ 2 ∇ 2 f Z (θ)) . ( 28
)
RW and Gibbs RW. For comparisons, we also run Random Walk-based samplers (RW) with Gaussian proposals: RW and Gibbs RW are defined respectively from algorithm 2 and algorithm 3, with PGdec and PGdual, γ 1 is the same step size as in ( 25) and [START_REF] Kent | Time-Reversible Diffusions[END_REF]. For the O-part (block j = 2), we choose C 1,2 := 2γ 2 I T ; for PGdec and PGdual, γ 2 is the same step size as in ( 26) and [START_REF] Roberts | Langevin Diffusions and Metropolis-Hastings Algorithms[END_REF].

I 1 = I 2 = 1 and µ 1,j (θ) = θ for j = 1, 2.
Twelve sampling strategies. The present section will compare twelve different samplers, constructed from the three drift strategies (PGdec, PGDual, RW), times two families (Metropolis-Hastings, Gibbs), times two choices for the covariance matrix C i,1 (I, O).

Step sizes. Different strategies are compared for the definition of the step sizes (γ 1 , γ 2 ). All of them consist in adapting the step sizes during the burn-in phase in such a way that the mean acceptance-rejection rate reaches approximately 0.25 (which is known to be optimal for some RW Metropolis-Hastings samplers, [START_REF] Gelman | Weak convergence and optimal scaling of random walk Metropolis algorithms[END_REF]). We observed that convergence occurs before 5e5 iterations, see Fig. 1[row 1, columns 2 and 3]; the proposed samples are all rejected and the chain does not move from its initial point during the first iterations, when the step size is too large. At the end of the burn-in phase, the step sizes are frozen and no longer adapted. For the Metropolis-Hastings samplers PGdec, PGdual, and RW, γ 1 is adapted and

γ 2 := γ 1 (case PGdec) of γ 2 := (λ R /λ O ) 2 γ 1 (cases PGdual, RW).
For the Gibbs samplers, the acceptancerejection steps are specific to each block R and O. A first consequence is that a move on one block can be accepted while the other one is not; this may yield larger step sizes (γ 1 , γ 2 ) which in turn favor larger moves of the chains. A second consequence is that we use the acceptance-rejection rate for the R-part (resp. for the O-part) to adapt γ 1 (resp. γ 2 ): the two step sizes have their own efficiency criterion.

C. Performance assessment criteria

The sampler performances are assessed and compared using three different criteria, see Fig. 1. Distance to the MAP. We compute the normalized distance R n -R MAP / R MAP where R MAP denotes the MAP estimator, (computed as in [START_REF] Pascal | Nonsmooth convex optimization to estimate the covid-19 reproduction number space-time evolution with robustness against low quality data[END_REF]): it is displayed vs. the iteration index n, in the burn-in phase (row 1) and after the burnin phase (row 2). This criterion quantifies the ability of the chains to perform a relevant exploration of the distribution: The chains have to visit the support of π Z , they show better ergodicity rates when they are able to rapidly escape from low density regions to move to higher density regions. Paths start from a non-informative initial point, considered as a point in a low density region. This criterion permits to quantify a relevant behavior of the Markov chain when (i) it drifts rapidly towards zero during the burn-in phase and (ii) it fluctuates in a large neighborhood of zero after the burn-in phase. Autocorrelation function (ACF). We then compute the ACF for each of the 2T components of the vector θ along the Markov path, for a lag from 1 to 1e5. On row 3, we report the mean value, over these 2T components, of the absolute value of the ACF versus the lag. This criterion quantifies a relevant behavior of the Markov chain when it converges rapidly to zero; it is indeed related to the effective sample size of a Markov chain (see [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]) and to the mixing properties of the chain (see e.g. [START_REF] Meyn | Markov Chains and Stochastic Stability[END_REF]Theorem 17.4.4 and Section 17.4.3.]). For example, a weaker ACF means that less iterations of the sampler are required to reach a given estimation accuracy. Gelman-Rubin (GR) statistic. Finally, we compute the GR statistic (see [START_REF] Gelman | Inference from Iterative Simulation Using Multiple Sequences[END_REF][START_REF] Brooks | General methods for monitoring convergence of iterative simulations[END_REF]) which quantifies a relevant behavior of the Markov chain when it converges rapidly to one. It measures how the sampler forgets its initial value and provides Fig. 1: Displayed: the normalized distance to R MAP versus the iteration index, during the burn-in phase (row 1) and after the burn-in phase (row 2); the ACF criterion versus the lag (row 3) and the GR criterion versus the iteration index (row 4); for the RW-based samplers (column 1), PGdec-based samplers (column 2) and PGdual-based samplers (column 3). For each Algo in { RW, PGdec, PGdual }: Algo I is in cyan, Algo O is in red, Gibbs Algo I is in blue and Gibbs Algo O is in black.

homogeneous approximations of the target distribution π Z after a given number of iterations. On row 4, we report this statistic versus the iteration index.

D. Performance comparisons.

Covariance matrices. Normalized Distance to MAP indices (Fig. 1[rows 1 & 2]) and GR indices (Fig. 1[row 4]) show that, for all algorithms RW, PGdec, and PGdual, the strategy O is more efficient than the strategy I, as corresponding indices decay more rapidly (to 0 and 1 respectively) for the formers than for the latters. A plausible explanation is that the smaller condition number of D O compared to D I is crucial to improve the mixing behavior of the Markov chains.

Metropolis-Hastings vs. Gibbs samplers. The evolution of the ACF criteria (Fig. 1[row 3]) shows that the Gibbs strategies are globally more efficient than the Metropolis-Hastings ones. Further, GR indices confirm better efficiency of the Markov chains obtained from Gibbs O strategies. Drift functions. The benefit of using functions µ miming optimization algorithms to drift the proposed points towards the higher probability regions of π Z is clearly quantified in Fig. 1 across all performance indices. During the burn-in phase, the normalized distance to MAP indices decay toward 0 significantly faster for all PGdualbased algorithms, compared to RW-based ones. Also, the PGdual O algorithms reach an expected plateau early in the burn-in phase, while this is barely the case at the end of the burn-in phase for RW O algorithms.

After the burn-in phase (second row), the samplers using optimization-based drift functions µ show better behavior after reaching the high density regions as they permit a broader and more rapid exploration of the support of the distribution around its maximum, with large amplitude moves from R MAP , and faster returns to R MAP . This is notably clearly visible for the Gibbs PGdual O strategy.

ACF indices decay more rapidly for the PGdual strategies than for the RW ones, irrespective of the choice of the covariance matrix or of the Metropolis-Hastings or Gibbs versions. Also, ACF indices for PGdual algorithms are less sensitive to the choice of the covariance matrices than the RW ones.

Finally, the GR statistic indices (row 4) clearly illustrates that Markov chains produced by PGdual samplers have better mixing properties, compared to others, showing hence sensitivity to the choice of the initial point. Optimal sampler. Combined together, these observations, globally consistent with those stemming from other countries or time periods, lead to the following generic comments.

While yielding essentially equivalent performance, across time periods and countries, the choice of the covariance matrices (algorithms O and I) significantly outperforms algorithms relying on Identity covariance matrices (not reported here as showing poor performances). This non trivial construction actually stemmed from the thorough and detailed mathematical analysis conducted in subsection III-C. Theorem 10 advocates to perform the augmentation of the (T -2) × T matrix D into a T × T invertible matrix D, by adding two rows which are orthogonal to the rows of D.

The observation that Gibbs samplers show better performances may stem from the fact that they benefit from larger values of the step sizes learned on the fly by the algorithms which favor larger jumps when proposing θ n+1/2 from θ n and imply better mixing properties.

Finally, the PGdual algorithms show homogeneous performances for example when varying the covariance matrix, and are thus less sensitive to parameter tuning, an important practical feature. Also, the PGdual algorithms show systematically better performances than the PGdec ones. This may result from the definition of µ PGdec which, because of the blocksplitting approach (see [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]), uses only partial information on π Z at each iteration.

As an overall conclusion, systematically observed across all studied time periods and countries, the Gibbs PGdual O algorithm, devised from the careful analysis of the theoretical properties of the distributions defined by A1-A2, is consistently found to be the most efficient strategy for a relevant assessment of the Covid19 pandemic time evolution.

VI. CREDIBILITY INTERVALS FOR R

Goal and set-up. The present section aims to illustrate the relevance of credibility interval-based estimations for the reproduction number R and for the outliers O from real Covid19 data. The statistical model for (R, O) is described in Section II-B. Because it is of greater interest for epidemiologists, the credibility intervals are translated into credibility intervals for the denoised counts Z (D) by their simple subtraction to the original counts Z, i.e., intuitively Z (D) = Z -O.

Credibility interval-based estimations are reported for Gibbs PGdual O sampling scheme only, as Section V established that it achieved the best performance amongst the twelve sampling schemes tested. Estimations are computed for a time period of five weeks (T = 35 days), which corresponds to a few typical pandemic time scales, of the order of 7 days, induced by the serial interval function Φ, cf. Section II. This period is set to a recent phase of the pandemic (Dec. 13th, 2021 to Jan. 17th, 2022). Estimations are reported for several countries, arbitrarily chosen as representatives of the pandemic across the world, but conclusions are valid for most countries. Computation of the credibility intervals. Credibility intervals are computed as follows: For each day t ∈ [T 1 , T 2 ] in the period of interest [T 1 , T 2 ], the chosen sampling scheme (Gibbs PGdual O) outputs 7.10 6 points of a Markov chain approximating the distribution of R t ; Quantiles of that distribution are estimated using the empirical cumulative distribution function ; For a chosen credibility level of 1α, the upper and lower limits of the credibility intervals are defined by the empirical 1α/2 and α/2 quantiles. For illustration purposes, 1α is set here to 95%. The same procedure is applied to produce credibility intervals for the outliers O t . Finally, credibility intervals for the estimated denoised new infection counts Z (D) are obtained by subtracting the credibility intervals for O t to the count Z t . Fig. 2 reports, top plots, the daily counts of new infections (black lines) to which are superimposed the 95% credibility interval-based estimations for the denoised counts, Z (D) , (red pipes). Fig. 2 further reports, bottom plots, the 95% credibility interval-based estimations for R. Relevance of credibility interval-based estimations. Fig. 2, together with the examination of equivalent plots for other countries, yields the following generic conclusions.

The estimated denoised counts show far smoother evolution along time compared to the raw counts, hence providing far more realistic assessments of the intensity of the pandemics. Notably, for most countries, the zero or low counts, associated with week-ends or non-working days, followed by days with over-evaluated counts by compensation, are smoothed out by the outlier estimation procedure, while the values of the counts for the regular (or non corrupted) days are left unchanged. This is the direct benefit of the nonlinear filtering procedure underlying the estimation formulation in [START_REF] Flahault | COVID-19 cacophony: is there any orchestra conductor?[END_REF], as opposed to traditional denoising procedures performed by classical linear filtering (such as moving average).

For the credibility intervals for R, their sizes range, depending on countries and time periods, from below 1% to above 10%. A careful examination shows that the credibility interval size is mostly driven by data quality: the credibility interval size increases when outliers are detected. Further, for a given country, the size of the credibility intervals varies only mildly along time over a five-week period. Changes in size are often associated with changes in the trends of the estimates of R, or with the occurrence of outliers. These credibility intervals provide hence a relevant assessment not only of the intensity of the pandemics, but also of the confidence that can be granted to this assessment, by providing epidemiologists with a range of likely values of R, rather than a single value. This permits to compare the evolution of the pandemics across several countries on a better scientifically grounded basis.

These credibility interval-based estimations permit a double analysis of the pandemic: They permit retrospectively to evaluate the impacts of sanitary measures on the pandemic evolution. Additionally, the smooth nature of the estimation of R (close to piecewise linear) performs an implicit short term forecast (or a nowcast) of the evolution of the pandemic intensity: For instance, for several countries (e.g., France, Mali, Brazil, Singapore,. . . ) the estimate of R is decreasing for the last 5 to 10 days of the studied period, predicting that daily new infections will reach a maximum of the current wave within the coming days and then will start to decrease.

These credibility interval estimates can be complemented with other estimates such as the Maximum, median or Mean a Posteriori (cf. e.g., [START_REF] Artigas | Credibility interval design for covid19 reproduction number from nonsmooth langevin-type monte carlo sampling[END_REF]63]).

Finally, let us emphasize that these estimates, denoised counts and credibility intervals, are obtained using a single and same set of hyperparameters λ R and λ O common to all countries.

VII. CONCLUSIONS

The proposed tools perform a relevant credibility intervalbased estimation of the Covid19 reproduction number and denoised new infection counts, by combining a statistical modeling of the time evolution of the pandemic with Monte Carlo Metropolis sampling strategies. Robustness against the low quality of the Covid19 is achieved by engineering the Bayesian model to impose sparsity in the changes of a smooth time evolution for the reproduction number and in the outlier occurrences, modeling data corruption. This is obtained at the price of the non-differentiability of their a posteriori distribution, thus precluding the use of the classical Metropolis Adjusted Langevin Algorithm to produce credibility intervals. This lead us to propose several Proximal-Gradient Monte Carlo algorithms tailored to the sampling of non-differentiable distributions, that also constitute valid sampling schemes for a much broader range of applications than that of the strict Covid19 pandemic monitoring, e.g., in image processing or more generally, in any Bayesian inverse problems with several nonsmooth priors.

Estimation performances were assessed and compared on real Covid19 data, using a set of well-selected indices quantifying the efficiency of these different sampling schemes.

Finally, it was shown for several countries and for a recent five-week period that the achieved credibility intervalbased estimations of both denoised new infection counts and reproduction number provide practitioners with an efficient and robust tool for the actual and practical monitoring of the Covid19 pandemic. Such estimates are updated on a daily basis on the authors's web-pages. Automated and data-driven estimations of the hyperparameters λ R and λ O are under current investigations.

In an effort toward reproducible research and open science, MATLAB codes implementing PGdec, PGdual for both Metropolis-Hastings and Gibbs versions are made publicly available at https://github.com/gfort-lab/OpSiMorE.

Second statement. Let {( R n , O n ), n ≥ 0} be a sequence in D Z and τ such that Z τ > 0. By (33), we have

C( R n , O n ) + t =τ ln(Z t !) ≥ I τ (D -1 R n , O n ) -Z τ ln I τ (D -1 R n , O n ).
The RHS tends to +∞ since

lim n I τ (D -1 R n , O n ) = 0 and Z τ > 0 by assumptions, hence, lim n C( R n , O n ) = +∞.
◮ Conclusion: a MAP exists. Assume that there exist t ⋆ < t ⋆⋆ in {1, . . . , T } such that Φ Z t⋆ > 0 and Φ Z t⋆⋆ > 0. Then ln π possesses at least one minimizer in D Z .

Proof. Consider the point (R 0 , O 0 ) given by R 0 := (1, . . . , 1) ⊤ and O

0 := (1, . . . , 1) ⊤ . Then, (R 0 , O 0 ) ∈ D Z . Set M 0 := -ln π(R 0 , O 0 ).
The goal of the proof below is to build a closed bounded set K in D Z such that outside K,ln π ≥ 1 + M 0 . This implies that (R 0 , O 0 ) ∈ K. Sinceln π is continuous on D Z , it reaches its minimum on the compact subset K, and this minimum is upper bounded by M 0 . Hence, this minimizer is also a global minimizer. Let us define K. We have

λ min R 2 ≤ D -1 R 2 = R 2 ≤ λ max R 2
where λ min (resp. λ max ) is the minimal (resp. maximal) eigenvalue of D -⊤ D -1 ; they are positive and finite. Consequently, by setting R

n := DR n , we have ( R n , O n ) ℓ → +∞ iff (R n , O n ) ℓ → +∞ for ℓ = 1, 2
since the norms are equivalent on R 2T . This property, the coercivity property (statement 1), and the equality [START_REF] Dwivedi | Logconcave sampling: Metropolis-Hastings algorithms are fast![END_REF]

imply that lim n -ln π(R n , O n ) = +∞ for any D Z -valued sequence {(R n , O n ), n ≥ 0} such that lim n R n + lim n O n = +∞. As a consequence, there exists C 1+M 0 such that (R, O) ∈ D Z , R + O > C 1+M0 =⇒ -ln π(R, O) ≥ 1 + M 0 . Similarly, there exists c 1+M0 > 0 such that (R, O) ∈ D Z , R t Φ Z t + O t < c 1+M 0 for some t s.t. Z t > 0 =⇒ -ln π(R, O) ≥ 1 + M 0 .
Consequently, we define

K := D Z ∩ {θ : R + O ≤ C 1+M0 } ∩ {θ : R t Φ Z t + O t ≥ c 1+M 0 for t s.t. Z t > 0}.

C. About the uniqueness of the MAP

We just proved that there exists a compact subset of the interior of D Z that contains a minimizer ofln π. Let θ ⋆ = (R ⋆ , O ⋆ ) be a minimizer.

• One or uncountably many. The functionln π is convex and finite on a convex set: hence, given a second minimizer θ ⋆⋆ , µθ ⋆ +(1-µ)θ ⋆⋆ is also a minimizer, whatever µ ∈ [0, 1].

• Same intensity, data fidelity term, and penalty term. Let f Z and g be defined by ( 2), (4). Following the same lines as in [START_REF] Pascal | Nonsmooth convex optimization to estimate the covid-19 reproduction number space-time evolution with robustness against low quality data[END_REF] where the strict convexity of the Kullback-Leibler term f Z is the key ingredient, it can be proved that I t (θ ⋆ ) = I t (θ ⋆⋆ ) for any t ∈ {1, . . . , T } and thus f Z (θ ⋆ ) = f Z (θ ⋆⋆ ); since ln π(θ ⋆ ) =ln π(θ ⋆⋆ ) since both points minimizeln π, we have g(θ ⋆ ) = g(θ ⋆⋆ ).

• Sign conditions. Set sign(a) = 1 when a > 0, sign(a) = -1 when a < 0 and sign(a) = 0 when a = 0. Set θ µ := θ ⋆ + µ(θ ⋆⋆θ ⋆ ); as proved above, f Z (θ µ ) = f Z (θ ⋆ ). We prove that if the sign conditions do not hold, for µ small enough g(θ µ ) < g(θ ⋆ ) which yields a contradiction since θ ⋆ is a minimizer. By [START_REF] Luu | Sampling from Non-smooth Distributions Through Langevin Diffusion[END_REF], we obtain For µ small enough, the RHS is negative when for some t, sign(O ⋆ t )sign(O ⋆⋆ t ) = -1 or when for some s, sign((DR ⋆ ) s )sign((DR ⋆⋆ ) s ) = -1. If such, g(θ µ ) < g(θ ⋆ ).

O ⋆ +µ(O ⋆⋆ -O ⋆ ) 1 -(1-µ) O ⋆ 1 -µ O ⋆⋆ 1 = µ T t=1 (
• Sufficient conditions for uniqueness. The proof is adapted from [START_REF] Ali | The Generalized Lasso Problem and Uniqueness[END_REF]Section 4]. Let θ ⋆⋆ = θ ⋆ + ω be another minimizer. Define

U := λ R D 0 (T -2)×T 0 T ×T λ O I T ∈ R (2T -2)×(2T ) .
The Fermat rule ([24, Theorem 16.2]) which characterizes optimality implies that zero is in the subdifferential ofln π at θ ⋆ : there exists γ(θ ⋆ ) ∈ R 2T -2 such that ∇f Z (θ ⋆ ) + U ⊤ γ(θ ⋆ ) = 0 2T ×1 where γ(θ ⋆ ) is the subgradient of the L 1norm in R 2T -2 evaluated at Uθ ⋆ ∈ R 2T -2 . Since ∇f Z (θ ⋆ ) depends on θ ⋆ through the I t 's which are constant on the set of the minimizers M (see above), and U ⊤ has full column rank, γ(θ ⋆ ) is the same whatever the minimizer θ ⋆ ; it is denoted by γ ⋆ . Set I := {j ∈ {1, . . . , 2T -2} : |γ ⋆ j | < 1}. Observe that any minimizer is in the kernel K 1 of the matrix U I which, by definition, collects the rows of U indexed by I; hence ω ∈ K 1 . In addition, ω is in the kernel K 2 of the T × (2T ) matrix [diag(Φ Z 1 , . . . Φ Z T ) I T ], since all the I t 's are constant on M. Therefore, if K 1 ∩ K 2 = {0}, the MAP is unique. D. Since π < ∞ on D then sup K π < ∞ by [B2]. For any measurable set A ⊆ K and any θ ∈ K, it holds

P ⋆ (θ, A) ≥ ι∈S ρ ι (θ) A q ι (θ, θ ′ )α ι (θ, θ ′ )dθ ′ ≥ A ι∈S ρ ι (θ)q ι (θ, θ ′ ) π(θ ′ ) ∧ ρ ι (θ ′ )q ι (θ ′ , θ) π(θ) π(θ ′ )dθ ′ ≥ inf K×K ι∈S ρ ι q ι sup K π A π(θ ′ )dθ ′ .
The RHS is positive by [B3] and this proves that K is 1-small and the chain is aperiodic.

• Unique invariant probability distribution. where ι = (i 1 , • • • , i J ) and µ is µ PGdec or µ PGdual . Therefore, since µ i,j (τ ) < ∞ for any τ ∈ D, we have q ι (θ, θ ′ ) ∧ q ι (θ ′ , θ) > 0 for any ι ∈ S and θ, θ ′ ∈ D.

• Both algorithms satisfy [B3]. For any compact set K of D, we have sup K µ i,j < ∞; in addition, µ is a continuous function on D (the function f is continuously differentiable and the proximal operator is continuous by [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Proposition 12.28]. Hence inf K×K q ι > 0 and [B3] holds.

• Positive Harris recurrence. From Claim 1, the PGdec Markov chain and the PGdual one are positive recurrent (they are φ-irreducible with an invariant distribution, and recurrent by [START_REF] Meyn | Markov Chains and Stochastic Stability[END_REF]Proposition 10.4.4]). Following the same lines as in [START_REF] Roberts | Harris recurrence of Metropolis-within-Gibbs and trans-dimensional Markov chains[END_REF]Theorem 8], we prove that the chain is Harris recurrent by showing that for any measurable set A such that A π(θ)dθ = 1 and any θ ∈ D, P θ (τ A < ∞) = 1 where τ A is the return-time to the set A ([54, Theorem 6(v)]); here P θ denotes the probability on the canonical space of the Markov chain with initial distribution the Dirac mass at θ and with kernel P ⋆ . Let A be a measurable subset of D such that A π(θ)dθ = 1. Let θ ∈ D. We write the kernel P ⋆ as follows P ⋆ (θ, A) = (1r(θ))M (θ, A) + r(θ)δ θ (A), where r(θ) := 1 -ι∈S ρ ι (θ) D q ι (θ, θ ′ )α ι (θ, θ ′ )dθ ′ , and M (θ, A) := (1r(θ)) -1 ι∈S ρ ι (θ) A q ι (θ, θ ′ )α ι (θ, θ ′ )dθ ′ .

Hence, P ⋆ (θ, •) is a mixture of two distributions: a Dirac mass at θ and M (θ, •). Since A c π(θ)dθ = 0 (here, A c := D \ A), then the Lebesgue measure of A c is 0. This implies that M (θ, A c ) = 0 and M (θ, A) = 1. It holds P θ (τ A = +∞) = E θ [½ X1 / ∈A P X1 (τ A = +∞)] = E θ [½ X1∈A c P X1 (τ A = +∞)] = r(θ) P θ (τ A = +∞); indeed, starting from θ, the chain can not reach A c when the kernel M (θ, •) is selected; and remains at θ when this kernel is not selected. Since r(θ) < 1 (otherwise the chain can not be φ-irreducible), we have P θ (τ A = +∞) = 0. This concludes the proof.

Example 11 (

 11 Example 2 and Example 9, to follow). A first strategy is to decompose g as in Example 2; Example 9 provides a possible augmentation of D. Another one is proposed by Theorem 10:

Proposition 13 .

 13 Assume A1, A2 and A3. Assume also that π is continuous on D. Then the sequence {θ n , n ≥ 0} given by algorithm 2 applied with µ = µ PGdec is a Markov chain, taking values in D and with unique invariant distribution π.In addition, for any initial point θ 0 ∈ D and any measurable functionh such that |h(θ)| π(θ)dθ < ∞, n ) = h(θ) π(dθ),with probability one. The same result holds when µ = µ PGdual . Proof. We have θ n+1 ∈ D since π(θ n+1/2 ) = 0 when θ n+1/2 / ∈ D. The other properties result from [52, Lemmas 1.1. and 1.2] and [53, Propositions 10.1.1 and 10.4.4 and Theorem 17.0.1]. The Harris recurrence property (required for the law of large numbers to hold for any θ 0 ) can be proved along the same lines as [54, Theorem 6(v)]. See the detailed proof in section IX of the Supplementary material.

Algorithm 3 : 5

 35 Gibbs Blockwise Proximal sampler Data: d j × d j positive definite matrices C i,j ; γ j > 0;a positive integer N max ; θ 0 ∈ D Result: A D-valued sequence {θ n , n ∈ {0, . . . , N max }} 1 Set ϑ = θ 0 ; 2 for n = 0, . . . , N max -1 do 3 for j = 1, . . . , J do 4 Sample i j ∈ {1, . . . , I j } with probability 1/I j ; Sample ξ n+1 j ∼ N dj (0 dj , C ij ,j ) ;

8

  Update ϑ = (θ n+1 1:j , θ n j+1:J )daily new infections and new death counts from the National Public Health Authorities of 200+ countries and territories of the world. Counts are updated on a daily basis since the early stage of the pandemic (Jan. 1st, 2020) until today. This repository thus provides researchers with a remarkable dataset to analyze the pandemic.

  ) where D := D I (see Example 9, the sampler is referred to as PGdual Invert (I)) or D := D O (see Example 11, the sampler is referred to as PGdual Ortho (O)); for the O-part, for both PGdual I and PGdual O,

  Covariance matrices C i,j . For the R-part (block j = 1), based on the comment in subsection III-C, we choose C 1,1 := 2γ 1 D I -1 D I -⊤ for PGdual I and Gibbs PGdual I; and C 1,1 := 2γ 1 D O -1 D O -⊤ for PGdual O and Gibbs PGdual O. The PGdec-based samplers and the RW-based samplers are run with the same covariance matrices. For PGdec, this yields C i,1 := 2γ 1 D I -1 D I -⊤ for PGdec I, Gibbs PGdec I and C i,1 := 2γ 1 D O -1 D O -⊤ for PGdec O, Gibbs PGdec O for the three cases i = 1, 2, 3. For RW, C 1,1 := 2γ 1 D I -1 D I -⊤ for RW I and Gibbs RW I; and C 1,1 := 2γ 1 D O -1 D O -⊤ for RW O and Gibbs RW O. For

Fig. 2 :

 2 Fig. 2: Gibbs PGdual O sampler-based estimation of the time evolution of R for a recent 5-week time period and for several countries. Top rows: Raw daily new infections counts Z (black) and estimates of the denoised counts Z (D) , obtained by subtracting the 95%-Credibility interval estimates of the outliers to the raw new infection counts Z ; Bottom rows: 95%-Credibility interval estimated of R t .

1 + µ d t=1 (

 1t=1 For A, B in R, and µ > 0, we have |A + µB| = |A| + µsign(A) sign(B)|B|½ A =0 + µ|B|½ A=0 + µo(1) where o(1) is a function satisfying lim µ→0 o(1) = 0. Hence, for τ = τ 1:d , τ ′ = τ ′ 1:d ∈ R d and µ ∈ (0, 1),τ + µ(τ ′τ ) 1 = (1µ) τ 1 + µ τ ′ sign(τ t )sign(τ ′ t ) -1) |τ ′ t |½ τt =0 + µ o[START_REF] Flahault | COVID-19 cacophony: is there any orchestra conductor?[END_REF].[START_REF] Luu | Sampling from Non-smooth Distributions Through Langevin Diffusion[END_REF]

Finally, [ 53 ,Claim 2 .

 532 Propositions 10.1.1. and 10.4.4] prove that π is the unique invariant distribution. Both PGdec and PGdual satisfy [B1,B2,B3]. The PGdec and PGdual chains are positive Harris-recurrent Markov chains: they satisfy a strong law of large numbers for any initial value in D.Proof. • Both algorithms satisfy [B1]. For both algorithms,ρ ι (θ) = 1/(I 1 • • • I J ) and q ι (θ, θ ′ ) is proportional to J j=1 exp -0.5(θ ′ jµ ij ,j (θ)) ⊤ C -1 ij ,j (θ ′ jµ ij ,j (θ))

•

  Strong Law of Large numbers. A positive Harris recurrent chain satisfies a strong law of large numbers whatever the initial value in D ([53, Theorem 17.0.1]).

  sign(O ⋆ t )sign(O ⋆⋆ t ) -1) |O ⋆⋆ t |½ O ⋆ t =0 + µ o(1).

	We T
	t=1	(sign(O ⋆ t ) sign(O ⋆⋆ t ) -1) |O ⋆⋆ t |½ O ⋆

have a similar expansion for

DR ⋆ + µ(DR ⋆⋆ -DR ⋆ ) 1 . Since g(θ ⋆ ) = g(θ ⋆⋆ ), this yields g(θ µ )g(θ ⋆ ) = µ T -2 t=1 |(DR ⋆⋆ ) t | . . . × (sign((DR ⋆ ) t ) sign((DR ⋆⋆ ) t ) -1) ½ (DR ⋆ )t =0

+ µ t =0 + µ o(1).

https://coronavirus.jhu.edu/

https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse covid 19 tim

Work partly supported by the Fondation Simone et Cino Del Duca, Institut de France.

Supplementary material VIII. PROOF OF PROPOSITION 1

A. Notations

Let D Z be the subset of (R + ) T × R T given by [START_REF] Obadia | The R0 package: A toolbox to estimate reproduction numbers for epidemic outbreaks[END_REF]. To shorten the notations, we will write I t (θ) instead of I t (θ, Z):

and π instead of π Z (see [START_REF] Flahault | COVID-19 cacophony: is there any orchestra conductor?[END_REF]) . Observe thatln π(θ) is equal to +∞ for θ / ∈ D Z and for θ ∈ D Z ln π(θ

for some normalizing constant C π . By convention 0 ln 0 = 0. Define the T × T invertible matrix D and compute its inverse:

Finally, define the criterion (30) for ( R, O) in the set

B. Existence of a MAP

We start with two lower bounds on the criterion C which will help us to study its behavior on some boundaries of

Proof. For any p > 0 and z ∈ N, p z exp(-p)/z! ∈ (0, 1) thus implying that z ln ppln(z!) ≤ 0. When p = z = 0, then z½ z>0 ln pp = 0 and ln(z!) = 0; hence we have

from which we obtain [START_REF] Roberts | Exponential convergence of Langevin distributions and their discrete approximations[END_REF] and [START_REF] Durmus | Analysis of Langevin Monte Carlo via Convex Optimization[END_REF].

◮ Behavior of C on some boundaries of D Z . Assume that there exist t ⋆ < t ⋆⋆ in {1, . . . , T } such that Φ Z t⋆ > 0 and

2) Let t ∈ {1, . . . , T } such that Z t > 0. For any sequence

For the discussions below, remember that this implies that lim

We distinguish two cases. First, either O n 1 tends to infinity or R n 3:T 1 tends to infinity. In the second case, these two norms are assumed bounded but R n 1:2 1 tends to infinity.

• Second case. Now consider the case when

By definition of D -1 , we have for any t ∈ {1, . . . , T },

by convention, the last term in the RHS is zero when t = 1, 2.

Under the assumptions of this second case,

for any t ∈ {3, . . . , T }. We prove that either (D -1 R n ) t⋆ tends to infinity, or (D -1 R n ) t⋆⋆ tends to infinity -which will imply by [START_REF] Durmus | Analysis of Langevin Monte Carlo via Convex Optimization[END_REF] that the criterion tends to infinity.

= +∞ (observe that this limit can not be -∞, as a consequence of the assumptions of Case 2, ( 35) and ( 34)). Apply [START_REF] Durmus | Efficient Bayesian Computation by Proximal Markov Chain Monte Carlo: When Langevin Meets Moreau[END_REF]

which implies that [START_REF] Durmus | Analysis of Langevin Monte Carlo via Convex Optimization[END_REF].

Since t ⋆ < t ⋆⋆ , this equality and [START_REF] Chatterji | Langevin Monte Carlo without smoothness[END_REF]

We then conclude, along the same lines as in Subcase 2A, that

IX. PROOFS OF SECTION III

A. Detailed proof of Theorem 10

Throughout the proof, we write A, U and Ā as a shorthand notation for A i,j , U i,j and Āi,j . Under the stated assumptions,

We first focus on the gradient step in ( 16) leading to:

Second, for any τ = τ 1:dj ∈ R dj ,

= τ 1:dj -ci,j prox γj gi,j (τ dj -ci,j+1:dj ) ,

since under the stated assumptions, we have

Therefore,

.

Now, let us apply Ā-1 ; by [START_REF] Atchadé | A Moreau-Yosida approximation scheme for a class of high-dimensional posterior distributions[END_REF], we have

Finally, since Ā-1 Ā = I dj , we have

and this concludes the proof of [START_REF] Pascal | Block-coordinate proximal algorithms for scale-free texture segmentation[END_REF].

When AA ⊤ = I ci,j , we have Ω i,j = A ⊤ A and A Ω i,j = A.

When UU ⊤ = I dj -ci,j , we have Ω i,j = UU ⊤ and I dj -Π i,j = U ⊤ U; this yields

leading to the final result of Theorem 10.

B. Detailed proof of Proposition 13

Let a finite set S of indices. For any θ ∈ D, let {ρ ι (θ), ι ∈ S} be a weight function: ι∈S ρ ι (θ) = 1 and ρ ι (θ) ≥ 0. Finally, for any ι ∈ S, let q ι (θ, θ ′ )dθ ′ be a Markov transition with respect to the Lebesgue measure on R d . PGdec and PGdual are special instances of algorithm 4. They correspond to the case S := {(i 1 , . . . , i J ), i j ∈ {1, . . . , I j }}; ρ ι (θ) = 1/(I 1 I 2 . . . I J ) for any θ; and to

q ι (y, x) q ι (x, y)

and θ n+1 = θ n otherwise.

inf K×K ι∈S ρ ι q ι > 0. Then the sequence {θ n , n ≥ 0} obtained by algorithm 4 is a Markov chain, taking values in D. It is φ-irreducible, strongly aperiodic and π is its unique invariant distribution.

Proof. • θ n ∈ D for any n. The proof is by induction on n. This property holds true for n = 0. Assume that θ n ∈ D. If θ n+1/2 / ∈ D, then π(θ n+1/2 ) = 0 and α i (θ n , θ n+1/2 ) = 0, so that θ n+1 = θ n and θ n+1 is in D. This concludes the induction.

• π is an invariant probability measure. Conditionally to ι and θ n , the distribution of θ n+1 is

Following the sames lines as in [59, Theorem 7.2] (details are omitted), the detailed balance condition with π can be established

hence π is invariant for P ⋆ .

• Irreducibility. By [B1], the chain is φ-irreducible (see [START_REF] Mengersen | Rates of convergence of the Hastings and Metropolis algorithms[END_REF]Lemma 1.1.]).

• Aperiodicity. Let us prove that the compact sets are 1-small and the chain is aperiodic; the proof is on the same lines as the proof of [