Christian Robert

George Casella
email: casella@stat.ufl.edu

Olivier Cappé

David Spiegelhalter

Alan Gelfand

Peter Green

Jun Liu

Sharon Mcgrayne

Peter Müller

Gareth Roberts

A History of Markov Chain Monte Carlo * -Subjective Recollections from Incomplete Data

In this note we attempt to trace the history and development of Markov chain Monte Carlo (MCMC) from its early inception in the late 1940's through its use today.

We see how the earlier stages of the Monte Carlo (MC, not MCMC) research have led to the algorithms currently in use. More importantly, we see how the development of this methodology has not only changed our solutions to problems, but has changed the way we think about problems.

Introduction

Markov Chain Monte Carlo (MCMC) methods have been around for almost as long as Monte Carlo techniques, even though their impact on Statistics has not been truly felt until the very early 1990s, except in the specialized fields of spatial Statistics and image analysis where those methods appeared earlier.

(The emergence of Markov based techniques in Physics

What can be reasonably seen as the first MCMC algorithm is what we now call the Metropolis algorithm, published by [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF]. It emanates from the same group of scientists who produced the Monte Carlo method, namely the research scientists of Los Alamos, mostly physicists working on mathematical physics and the atomic bomb.2 MCMC algorithms therefore date back to the same time as the development of regular (MC only) Monte Carlo methods, which are usually traced to Ulam and von Neumann in the late 1940s. Stanislaw Ulam associates the original idea with an intractable combinatorial computation he attempted in 1946 (calculating the probability of winning at the card game "solitaire"). This idea was enthusiastically adopted by John von Neumann for implementation with direct applications to neutron diffusion, the name "Monte Carlo" being suggested by Nicholas Metropolis. [START_REF] Eckhardt | Stan ulam, john von neumann, and the monte carlo method[END_REF] describes these early Monte Carlo developments, and Hitchcock 2003 gives a brief history of the Metropolis algorithm.)

These occurrences very closely coincide with the appearance of the very first computer, the ENIAC, which came to life in February 1946, after three years of construction. The Monte Carlo method was set up by von Neumann, who was using it on thermonuclear and fission problems as early as 1947. At the same time, that is, 1947, Ulam and von Neumann invented inversion and accept-reject techniques (also recounted in [START_REF] Eckhardt | Stan ulam, john von neumann, and the monte carlo method[END_REF] to simulate from non-uniform distributions. Without computers, a rudimentary version invented by Fermi in the 1930s did not get any recognition [START_REF] Metropolis | The beginning of the Monte Carlo method[END_REF]. Note also that, as early as 1949, a symposium on Monte Carlo was supported by Rand, NBS and the Oak Ridge laboratory and that [START_REF] Metropolis | The Monte Carlo method[END_REF] published the very first paper about the Monte Carlo method.

2.1 The [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF]

paper

The first MCMC algorithm is associated with a second computer, called MANIAC(!), built 3 in Los Alamos under the direction of Metropolis in early 1952. Both a physicist and a mathematician, Nicolas Metropolis, who died in Los Alamos in 1999, came to this place in April 1943 . The other members of the team also came to Los Alamos during those years, with Edward Teller being the most controversial character of the group. As early as 1942, he was one of the first scientists to work on the Manhattan Project that led to the production of the A bomb. Almost as early, he became obsessed with the hydrogen (H) bomb, which he eventually managed to design with Stanislaw Ulam using the better computer facilities in the early 1950s. 4Published in June 1953 in the Journal of Chemical Physics, the primary focus of [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF] is the computation of integrals of the form I = F (p, q) exp{-E(p, q)/kT }dpdq exp{-E(p, q)/kT }dpdq , with the energy E being defined as

E(p, q) = 1 2 N i=1 N j=1 j =i V (d ij ),
where N is the number of particles, V a potential function and d ij the distance between particles i and j. The Boltzmann distribution exp{-E(p, q)/kT } is parameterized by the temperature T , k being the Boltzmann constant, with a normalization factor Z(T ) = exp{-E(p, q)/kT }dpdq that is not available in closed form. Since p and q are 2N-dimensional vectors, numerical integration is impossible. Given the large dimension of the problem, even standard Monte Carlo techniques fail to correctly approximate I, since exp{-E(p, q)/kT } is very small for most realizations of the random configurations of the particle system (uniformly in the 2N or 4N square). In order to improve the efficiency of the Monte Carlo method, [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF] propose a random walk modification of the N particles. That is, for each particle i

(1 ≤ i ≤ N), values x ′ i = x i + αξ 1i and y ′ i = y i + αξ 2i
are proposed, where both ξ 1i and ξ 2i are uniform U(-1, 1). The energy difference ∆E between the new configuration and the previous one is then computed and the new configuration is accepted with probability min {1, exp(-∆E/kT )} ,

and otherwise the previous configuration is replicated (in the sense that it will count one more time in the final average of the F (p t , p t )'s over the τ moves of the random walk, 1 ≤ t ≤ τ )).

Note that [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF] move one particle at a time, rather than moving all of them together, which makes the initial algorithm appear as a primitive kind of Gibbs sampler (!).

The authors of [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF] demonstrate the validity of the algorithm by first establishing irreducibility (that they call ergodicity) and second proving ergodicity, that is, convergence to the stationary distribution. The second part is obtained via a discretization of the space: They first note that the proposal move is reversible, then establish that exp{-E/kT } is invariant. The result is therefore proven in its full generality (modulo the discretization). The remainder of the paper is concerned with the specific problem of the rigid-sphere collision model. The number of iterations of the Metropolis algorithm seems to be limited: 16 steps for burn-in and 48 to 64 subsequent iterations (which still required four to five hours on the Los Alamos MANIAC).

An interesting variation of (1) is the Simulated Annealing algorithm, developed by [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF], who connected optimization with annealing, the cooling of a metal.

Their variation is to allow T of (1) to change as the algorithm runs, according to a "cooling schedule", and the Simulated Annealing algorithm can be shown to find the global maximum with probability 1, although the analysis is quite complex due to the fact that, with varying T , the algorithm is no longer a time-homogeneous Markov chain.

The Hastings (1970) paper

The Metropolis algorithm was later generalized by [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their application[END_REF] and [START_REF] Peskun | Optimum Monte Carlo sampling using Markov chains[END_REF][START_REF] Peskun | Guidelines for chosing the transition matrix in Monte Carlo methods using Markov chains[END_REF] as a statistical simulation tool that could overcome the curse of dimensionality met by regular Monte Carlo methods (already emphasized in [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF]).5

In his Biometrika paper, 6 Hastings (1970) also defines his methodology on finite and reversible Markov chains, treating the continuous case by using a discretization analogy.

The generic probability of acceptance for a move from state i to state j is

α ij = s ij 1 + π i π j q ij q ji
, where s ij is a symmetric function. This generic form of probability encompasses the forms of both [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF] and [START_REF] Barker | Monte Carlo calculations of the radial distribution functions for a proton electron plasma[END_REF]. At this stage, Peskun's ordering is not yet discovered and Hastings thus mentions that little is known about the relative merits of those two choices (even though) Metropolis's method may be preferable. He also warns against high rejection rates as indicative of a poor choice of transition matrix, but does not mention the opposite pitfall of low rejection rates, associated with a slow exploration of the target.

The examples given in the paper are a Poisson target with a ±1 random walk proposal, a normal target with a uniform random walk proposal mixed with its reflection (i.e. centered at -X(t) rather than X(t)), and then a multivariate target where Hastings introduces a Gibbs sampling strategy, updating one component at a time and defining the composed transition as satisfying the stationary condition because each component does leave the target invariant! [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their application[END_REF] actually refers to [START_REF] Erhman | Computation of order parameters in an ising lattice by the monte carlo method[END_REF] as a preliminary if specific instance of this sampler. More precisely, this is Metropolis-within-Gibbs except for the name. It is quite amazing that this first introduction of the Gibbs sampler has been completely overlooked, even though the proof of convergence is completely general, based on a composition argument as in [START_REF] Tierney | Markov chains for exploring posterior distributions (with discussion)[END_REF]! The remainder of the paper deals with (a) an importance sampling version of MCMC, (b) general remarks about assessment of the error, and (c) an application to random orthogonal matrices (with yet again an example of Gibbs sampling).

Three years later, again in Biometrika, [START_REF] Peskun | Optimum Monte Carlo sampling using Markov chains[END_REF] published a comparison of Metropolis' and Barker's forms of acceptance probabilities and shows (in a discrete setup) that the optimal choice (in terms of the asymptotic variance of any empirical average) is that of Metropolis. The proof is a direct consequence of a result by [START_REF] Kemeny | Finite Markov Chains[END_REF] on the asymptotic variance. Peskun also establishes that this asymptotic variance can improve upon the iid case if and only if the eigenvalues of P -A are all negative, when A is the transition matrix corresponding to the iid simulation and P the transition matrix corresponding to the Metropolis algorithm, but he concludes that the trace of P -A is always positive.

Seeds of the Revolution

A number of earlier pioneers had brought forward the seeds of Gibbs sampling; in particular, Hammersley and Clifford had produced a constructive argument in 1970 to recover a joint distribution from its conditionals, a result later called the Hammersley-Clifford theorem by [START_REF] Besag | Spatial interaction and the statistical analysis of lattice systems[END_REF][START_REF] Besag | On the statistical analysis of dirty pictures[END_REF]. Besides [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their application[END_REF] and [START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF], already mentioned, other papers that contained the germs of Gibbs sampling are [START_REF] Besag | Generalized Monte Carlo significance tests[END_REF], [START_REF] Broniatowski | Reconnaissance de mélanges de densités par un algorithme d'apprentissage probabiliste[END_REF], [START_REF] Qian | Parameter estimation for hidden Gibbs chains[END_REF], and [START_REF] Tanner | The calculation of posterior distributions by data augmentation[END_REF].

Besag's Early Work and the Fundamental (Missing) Theorem

In the early 1970's, Hammersley, Clifford, and Besag were working on the specification of joint distributions from conditional distributions and on necessary and sufficient conditions for the conditional distributions to be compatible with a joint distribution. What is now known as the Hammersley-Clifford theorem states that a joint distribution for a vector associated with a dependence graph (edge meaning dependence and absence of edge conditional independence) must be represented as a product of functions over the cliques of the graphs, that is, of functions depending only on the components indexed by the labels in the clique (which is a subset of the nodes of the graphs such that every node is connected by an edge to every other node in the subset). See [START_REF] Cressie | Spatial Statistics[END_REF] or [START_REF] Lauritzen | Graphical Models[END_REF] for detailed treatments.

From an historical point of view, [START_REF] Hammersley | Discussion of Mr Besag's paper[END_REF] explains why the Hammersley-Clifford theorem was never published as such, but only through [START_REF] Besag | Spatial interaction and the statistical analysis of lattice systems[END_REF]. The reason is that Clifford and Hammersley were dissatisfied with the positivity constraint: The joint density could be recovered from the full conditionals only when the support of the joint was made of the product of the supports of the full conditionals (with obvious counter-examples, as in [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]. While they strived to make the theorem independent of any positivity condition, their graduate student published [START_REF] Moussouris | Gibbs and Markov random systems with constraints[END_REF], a counter-example that put a full stop to their endeavors.

While Julian Besag can certainly be credited to some extent of the (re-)discovery of the Gibbs sampler (as in [START_REF] Besag | Spatial interaction and the statistical analysis of lattice systems[END_REF], [START_REF] Besag | Statistical analysis of non-lattice data[END_REF] has the curious and anticlimactic following comment:

The simulation procedure is to consider the sites cyclically and, at each stage, to amend or leave unaltered the particular site value in question, according to a probability distribution whose elements depend upon the current value at neighboring sites (...) However, the technique is unlikely to be particularly helpful in many other than binary situations and the Markov chain itself has no practical interpretation. So, while stating the basic version of the Gibbs sampler on a graph with discrete variables, Besag dismisses it as unpractical.

On the other hand, Hammersley, together with Handscomb, wrote a textbook on Monte Carlo methods, (the first?) [START_REF] Hammersley | Monte Carlo Methods[END_REF]. There they cover such topics as "Crude Monte Carlo" (which is (3)); importance sampling; control variates; and "Conditional Monte Carlo", which looks surprisingly like a missing-data completion approach. Of course, they do not cover the Hammersley-Clifford theorem but, in contrast to [START_REF] Besag | Spatial interaction and the statistical analysis of lattice systems[END_REF], they state in the Preface We are convinced nevertheless that Monte Carlo methods will one day reach an impressive maturity.

Well said!

EM and its Simulated Versions as Precursors

Besides a possible difficult computation in the E-step, problems with the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm (with discussion)[END_REF] do occur in the case of multimodal likelihoods. The increase of the likelihood function at each step of the algorithm ensures its convergence to the maximum likelihood estimator in the case of unimodal likelihoods but it implies a dependence on initial conditions for multimodal likelihoods. Several proposals can be found in the literature to overcome this problem, one of which we now describe because of its connection with Gibbs sampling. [START_REF] Broniatowski | Reconnaissance de mélanges de densités par un algorithme d'apprentissage probabiliste[END_REF] and Celeux andDiebolt (1985, 1992) have tried to overcome the dependence of EM methods on the starting value by replacing the E step with a simulation step, the missing data z being generated conditionally on the observation x and on the current value of the parameter θ m . The maximization in the M step is then done on the (simulated) complete-data log-likelihood, H(x, z m |θ). The appeal of this approach is that it allows for a more systematic exploration of the likelihood surface by partially avoiding the fatal attraction of the closest mode. Unfortunately, the theoretical convergence results for these methods are limited. [START_REF] Celeux | Une version de type recuit simulé de l'algorithme EM[END_REF] have, however, solved the convergence problem of SEM by devising a hybrid version called SAEM (for Simulated Annealing EM), where the amount of randomness in the simulations decreases with the iterations, ending up with an EM algorithm. This version actually relates to simulated annealing methods.

Gibbs, and Beyond

Although somewhat removed from statistical inference in the classical sense and based on earlier techniques used in Statistical Physics, the landmark paper by [START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF] brought Gibbs sampling into the arena of statistical application. This paper is also responsible for the name Gibbs sampling, because it implemented this method for the Bayesian study of Gibbs random fields which, in turn, derive their name from the physicist Josiah Willard Gibbs . This original implementation of the Gibbs sampler was applied to a discrete image processing problem and did not involve completion. But this was one more spark that led to the explosion, as it had a clear influence on Green, Smith, Spiegelhalter and others.

The extent to which Gibbs sampling and Metropolis algorithms were in use within the image analysis and point process communities is actually quite large, as illustrated in [START_REF] Ripley | Stochastic Simulation[END_REF] where Section §4.7 is entitled "Metropolis' method and random fields" and describes the implementation and the validation of the Metropolis algorithm in a finite setting with an application to Markov random fields and the corresponding issue of bypassing the normalizing constant. [START_REF] Besag | Bayesian image restoration, with two applications in spatial statistics (with discussion)[END_REF] is another striking example of the activity in the spatial statistics community at the end of the 1980's (the paper was submitted in 1989).

The Revolution

The gap of more than 30 years between [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF] and Gelfand and Smith (1990) can still be partially attributed to the lack of appropriate computing power, as most of the examples now processed by MCMC algorithms could not have been treated previously, even though the hundreds of dimensions processed in [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF] were quite formidable. However, by the mid-1980s, the pieces were all in place. After Peskun, MCMC in the statistical world was dormant for about 10 years, and then several papers appeared that highlighted its usefulness in specific settings like pattern recognition, image analysis or spatial statistics (see, for example, [START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF][START_REF] Tanner | The calculation of posterior distributions by data augmentation[END_REF][START_REF] Besag | Towards Bayesian image analysis[END_REF]). In particular, [START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF] building on [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF], [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their application[END_REF] , and [START_REF] Peskun | Optimum Monte Carlo sampling using Markov chains[END_REF], influenced Gelfand and Smith (1990) to write a paper that is the genuine starting point for an intensive use of MCMC methods by the (mainstream) statistical community. It sparked new interest in Bayesian methods, statistical computing, algorithms, and stochastic processes through the use of computing algorithms such as the Gibbs sampler and the Metropolis-Hastings algorithm.

(See Casella and George 1992 for an elementary introduction to the Gibbs sampler 7 .) Interestingly, the earlier paper by [START_REF] Tanner | The calculation of posterior distributions by data augmentation[END_REF] had essentially the same ingredients as Gelfand and Smith (1990), namely the fact that simulating from the conditional distributions is sufficient to simulate (in the limiting sense) from the joint. This paper was considered important enough to be a discussion paper in the Journal of the American Statistical Association, but its impact was somehow limited, compared with the one of Gelfand and Smith (1990). There are several reasons for this; one being that the method seemed to only apply to missing data problems (hence the name data augmentation), and another is that the authors were more focused on approximating the posterior distribution. They suggested a (Markov chain) Monte Carlo approximation to the target π(θ|x) at each iteration 7 On a humorous note, the original Technical Report of this paper was called Gibbs for Kids, which was changed because a referee did not appreciate the humor. However, our colleague Dan Gianola, an Animal Breeder at Wisconsin, liked the title. In using Gibbs sampling in his work, he gave a presentation in 1993 at the 44th Annual Meeting of the European Association for Animal Production, Arhus, Denmark. The title: Gibbs for Pigs. Looking at these meetings, we can see the paths that Gibbs sampling would lead us down. In the next two sections we will summarize some of the advances from the early to mid 1990s.

Advances in MCMC Theory

Perhaps the most influential MCMC theory paper of the 1990s is [START_REF] Tierney | Markov chains for exploring posterior distributions (with discussion)[END_REF], who carefully laid out all of the assumptions needed to analyze the Markov chains and then developed their properties, in particular, convergence of ergodic averages and central limit theorems.

In one of the discussions of that paper, [START_REF] Chan | Discussion of "Markov chains for exploring posterior distribution[END_REF] were able to relax a condition on Tierney's Central Limit Theorem, and this new condition plays an important role in research today (see Section 5.4). A pair of very influential, and innovative, papers is the work of [START_REF] Liu | Covariance structure of the Gibbs sampler with applications to the comparisons of estimators and sampling schemes[END_REF]Liu et al. ( , 1995)), who very carefully analyzed the covariance structure of Gibbs sampling, and were able to formally establish the validity of Rao-Blackwellization in Gibbs sampling. Gelfand and Smith (1990) had used Rao-Blackwellization, but it was not justified at that time, as the original theorem was only applicable to iid sampling, which is not the case in MCMC. Other early theoretical developments include the Duality Theorem of [START_REF] Diebolt | Estimation of finite mixture distributions by Bayesian sampling[END_REF], who showed that in the two-stage Gibbs sampler (which is equivalent to the Data Augmentation algorithm of [START_REF] Tanner | The calculation of posterior distributions by data augmentation[END_REF], convergence properties of one chain can be transferred to other chains, a fact also found in [START_REF] Liu | Covariance structure of the Gibbs sampler with applications to the comparisons of estimators and sampling schemes[END_REF]Liu et al. ( , 1995)). This turns out to be particularly important in mixture models, where it is typical that one part of the Gibbs chain is discrete and finite, and the other is continuous.

The convergence properties of the finite chain carry over to the continuous chain.

Another paper must be singled out, namely [START_REF] Mengersen | Rates of convergence of the Hastings and Metropolis algorithms[END_REF], for setting the tone for the study of the speed of convergence of MCMC algorithms to the target distribution.

Subsequent works in this area by Richard Tweedie, Gareth Roberts, Jeff Rosenthal and coauthors are too numerous to be mentioned here, even though the paper by Roberts et al.

(1997) must be cited for setting explicit targets on the acceptance rate of the random walk Metropolis-Hastings algorithm, as well as [START_REF] Roberts | Convergence of slice sampler Markov chains[END_REF] for getting an upper bound on the number of iterations (523) needed to approximate the target up to 1% by a slice sampler. The untimely death of Richard Tweedie in 2001 alas had a major impact on the book about MCMC convergence he was contemplating with Gareth Roberts.

One pitfall arising from the widespread use of Gibbs sampling was the tendency to spec-ify models only through their conditional distributions, almost always without referring to the positivity conditions in Section 3. Unfortunately, it is possible to specify a perfectly legitimate-looking set of conditionals that do not correspond to any joint distribution, and the resulting Gibbs chain cannot converge. [START_REF] Hobert | The effect of improper priors on Gibbs sampling in hierarchical linear models[END_REF] were able to document the conditions needed for a convergent Gibbs chain, and alerted the Gibbs community to this problem (which only arises if improper priors are used, but this is a frequent occurrence).

Much other work followed, and continues to grow today. [START_REF] Geyer | Annealing Markov chain Monte Carlo with applications to ancestral inference[END_REF] describe how to put a "ladder" of chains together to have both "hot" and "cold" exploration, followed by Neal's 1996 introduction of tempering; [START_REF] Athreya | On the convergence of the Markov chain simulation method[END_REF] gave more easily verifiable conditions for convergence; [START_REF] Meng | Seeking efficient data augmentation schemes via conditional and marginal augmentation[END_REF] and [START_REF] Liu | Parameter expansion for data augmentation[END_REF] developed the theory of parameter expansion in the Data Augmentation algorithm, leading to construction of chains with faster convergence, and to the work of [START_REF] Hobert | A theoretical comparison of the data augmentation, marginal augmentation and px-da algorithms[END_REF], who give precise constructions and theorems to show how parameter expansion can uniformly improve over the original chain.

Advances in MCMC Applications

The real reason for the explosion of MCMC methods was the fact that an enormous number of problems that were deemed to be computational nightmares now cracked open like eggs.

As an example, consider this very simple random effects model from Gelfand and Smith (1990). Observe

Y ij = θ i + ε ij , i = 1, . . . , K, j = 1, . . . , J, (2)
where

θ i ∼ N(µ, σ 2 θ ) ε ij ∼ N(0, σ 2 ε ), independent of θ i
Estimation of the variance components can be difficult for a frequentist (REML is typically preferred) but it indeed was a nightmare for a Bayesian, as the integrals were intractable.

However, with the usual priors on µ, σ 2 θ , and σ 2 ε , the full conditionals are trivial to sample from and the problem is easily solved via Gibbs sampling. Moreover, we can increase the number of variance components and the Gibbs solution remains easy to implement. During the early 1990s, researchers found that Gibbs, or Metropolis-Hastings, algorithms would crack almost any problem that they looked at, and there was a veritable flood of papers applying MCMC to previously intractable models, and getting good solutions. For example, building on (2), it was quickly realized that Gibbs sampling was an easy route to getting estimates in the linear mixed models [START_REF] Wang | Marginal inferences about variancecomponents in a mixed linear model using gibbs sampling[END_REF][START_REF] Wang | Bayesian analysis of mixed limear models via Gibbs sampling with an application to litter size in Iberian pigs[END_REF], and even generalized linear mixed models [START_REF] Zeger | Generalized linear models with random effects; a Gibbs sampling approach[END_REF]. Demarginalization (the introduction of latent variables) arguments made it possible to analyze probit models using a latent variable approach in a linear mixed model [START_REF] Albert | Bayesian analysis of binary and polychotomous response data[END_REF], and demarginalization was also a route to estimation in mixture models with Gibbs sampling (see, for example, [START_REF] Robert | Inference in mixture models[END_REF].

It progressively dawned on the community that latent variables could be artificially introduced to run the Gibbs sampler in about every situation, as eventually published in [START_REF] Damien | Gibbs sampling for Bayesian non-conjugate and hierarchical models by using auxiliary variables[END_REF], the main example being the slice sampler [START_REF] Neal | Slice sampling (with discussion)[END_REF]. A (very incomplete) list of some other applications include changepoint analysis [START_REF] Carlin | Hierarchical Bayesian analysis of change point problems[END_REF][START_REF] Stephens | Bayesian retrospective multiple-changepoint identification[END_REF]; Genomics [START_REF] Lawrence | Detecting subtle sequence signals -a Gibbs sampling strategy for multiple alignment[END_REF][START_REF] Stephens | Bayesian inference in multipoint gene mapping[END_REF][START_REF] Churchill | Accurate restoration of DNA sequences (with discussion)[END_REF]; capturerecapture (George andRobert 1992, Dupuis 1995); variable selection in regression [START_REF] George | Variable Selection Via Gibbbs Sampling[END_REF]; spatial statistics [START_REF] Raftery | Stopping the Gibbs sampler, the use of morphology, and other issues in spatial statistics[END_REF], and longitudinal studies [START_REF] Lange | Hierarchal bayes models for the progression of hiv infection using longitudinal cd4 t-cell numbers[END_REF]).

Many of these applications were advanced though other developments such as the Adaptive Rejection Sampling of [START_REF] Gilks | Derivative-free adaptive rejection sampling for Gibbs sampling[END_REF], [START_REF] Gilks | Adaptive rejection Metropolis sampling within Gibbs sampling[END_REF], and the simulated tempering approaches of [START_REF] Geyer | Annealing Markov chain Monte Carlo with applications to ancestral inference[END_REF] or [START_REF] Neal | Sampling from multimodal distributions using tempered transitions[END_REF].

After the Revolution

After the revolution comes the "second" revolution, but now we have a more mature field.

The revolution has slowed, and the problems are being solved in, perhaps, deeper and more sophisticated ways (even though Gibbs sampling also offers to the amateur the possibility to handle Bayesian analysis in complex models at little cost, as exhibited by the widespread use of BUGS). But, as before, the methodology continues to expand the set of problems that statisticians can provide meaningful solutions, and thus continues to further the impact of Statistics.

A Brief Glimpse at Particle Systems

The realization of the possibilities of iterating importance sampling is not new: in fact, it is about as old as Monte Carlo methods themselves! It can be found in the molecular simulation literature of the 50's, as in [START_REF] Hammersley | Poor man's Monte Carlo[END_REF], [START_REF] Rosenbluth | Monte Carlo calculation of the average extension of molecular chains[END_REF] and [START_REF] Marshall | The use of multi-stage sampling schemes in Monte Carlo computations[END_REF]. Hammersley and colleagues proposed such a method to simulate a self-avoiding random walk [START_REF] Madras | The Self-Avoiding Random Walk. Probability and its Applications[END_REF] on a grid, due to huge inefficiency in regular importance sampling and rejection techniques. Although this early implementation occurred in particle physics, the use of the term "particle" only dates back to [START_REF] Kitagawa | Monte Carlo filter and smoother for non-Gaussian non-linear state space models[END_REF], while [START_REF] Carpenter | Building robust simulation-based filters for evolving datasets[END_REF] coined the term "particle filter". In signal processing, early occurrences of a "particle filter" can be traced back to [START_REF] Handschin | Monte Carlo techniques to estimate the conditional expectation in multi-stage nonlinear filtering[END_REF].

More in connection with our theme, the landmark paper of [START_REF] Gordon | A novel approach to non-linear/non-Gaussian Bayesian state estimation[END_REF] introduced the bootstrap filter which, while formally connected with importance sampling, involves past simulations and possible MCMC steps [START_REF] Gilks | Following a moving target-Monte Carlo inference for Bayesian models[END_REF]. As described in the volume edited by [START_REF] Doucet | Sequential Monte Carlo Methods in Practice[END_REF], particle filters are simulation methods adapted to sequential settings where data are collected progressively in time as in radar detection, telecommunication correction or financial volatility estimation. Taking advantage of state-space representations of those dynamic models, particle filter methods produce Monte Carlo approximations to the posterior distributions by propagating simulated samples whose weights are actualized against the in-coming observations. Since the importance weights have a tendency to degenerate (that is, all weights but one are close to zero), additional MCMC steps can be introduced at times to recover the variety and representativeness of the sample.

Modern connections with MCMC in the construction of the proposal kernel are to be found, for instance, in [START_REF] Doucet | On sequential Monte Carlo sampling methods for Bayesian filtering[END_REF] and in Del [START_REF] Del Moral | The sequential Monte Carlo samplers[END_REF]. In parallel, sequential imputation was developed in [START_REF] Kong | Sequential imputations and Bayesian missing data problems[END_REF], while [START_REF] Liu | Blind deconvolution via sequential imputations[END_REF] first formally pointed out the importance of resampling in sequential Monte Carlo, a term coined by them.

The recent literature on the topic more closely bridges the gap between sequential Monte Carlo and MCMC methods by making adaptive MCMC a possibility (see, for example, Andrieu et al. 2004 or Roberts and[START_REF] Roberts | Coupling and ergodicity of adaptive mcmc[END_REF].

Perfect sampling

Introduced in the seminal paper of [START_REF] Propp | Exact sampling with coupled Markov chains and applications to statistical mechanics[END_REF], perfect sampling, namely the ability to use MCMC methods to produce an exact (or perfect) simulation from the target, maintains a unique place in the history of MCMC methods. Although this exciting discovery led to an outburst of papers, in particular in the large body of work of Møller and coauthors, including the book by [START_REF] Møller | Statistical Inference and Simulation for Spatial Point Processes[END_REF], as well as many reviews and introductory materials, like [START_REF] Casella | Explaining the perfect sampler[END_REF], [START_REF] Fismen | Exact simulation using Markov chains[END_REF][START_REF] Dimakos | A guide to exact simulation[END_REF], the excitement quickly dried out. The major reason for this ephemeral lifespan is that the construction of perfect samplers is most often close to impossible or impractical [START_REF] Foss | Perfect simulation and backward coupling[END_REF], despite some advances in the implementation (Fill 1998a,b).

There is, however, ongoing activity in the area of point processes and stochastic geometry, much from the work of Møller and Kendall. In particular, [START_REF] Kendall | Perfect simulation using dominating processes on ordered spaces, with application to locally stable point processes[END_REF] developed an alternative to the Coupling From The Past (CFPT) algorithm of [START_REF] Propp | Exact sampling with coupled Markov chains and applications to statistical mechanics[END_REF], called horizontal CFTP, which mainly applies to point processes and is based on continuous time birth-and-death processes. See also [START_REF] Fernández | Perfect simulation for interacting point processes, loss networks and Ising models[END_REF] for another horizontal CFTP algorithm for point processes. [START_REF] Berthelsen | Likelihood and non-parametric Bayesian MCMC inference for spatial point processes based on perfect simulation and path sampling[END_REF] exhibited a use of these algorithms for nonparametric Bayesian inference on point processes.

Reversible jump and variable dimensions

From many viewpoints, the invention of the reversible jump algorithm in [START_REF] Green | Reversible jump MCMC computation and Bayesian model determination[END_REF] can be seen as the second MCMC revolution: the formalization of a Markov chain that moves across models and parameters spaces allowed for the Bayesian processing of a wide variety of new models and contributed to the success of Bayesian model choice and subsequently to its adoption in other fields. There exist earlier alternative Monte Carlo solutions like [START_REF] Gelfand | Bayesian model choice: asymptotics and exact calculations[END_REF] and [START_REF] Carlin | Bayesian model choice through Markov chain Monte Carlo[END_REF], the later being very close in spirit to reversible jump MCMC (as shown by the completion scheme of [START_REF] Brooks | Efficient construction of reversible jump Markov chain Monte Carlo proposal distributions (with discussion)[END_REF], but the definition of a proper balance condition on cross-model Markov kernels in [START_REF] Green | Reversible jump MCMC computation and Bayesian model determination[END_REF] gives a generic setup for exploring variable dimension spaces, even when the number of models under comparison is infinite. The impact of this new idea was clearly perceived when looking at the First European Conference on Highly Structured Stochastic Systems that took place in Rebild, Denmark, the next year, organized by Stephen Lauritzen and Jesper Møller: a large majority of the talks were aimed at direct implementations of RJMCMC to various inference problems. The application of RJMCMC to mixture order estimation in the discussion paper of [START_REF] Richardson | On Bayesian analysis of mixtures with an unknown number of components (with discussion)[END_REF] ensured further dissemination of the technique. More recently, [START_REF] Stephens | Bayesian analysis of mixture models with an unknown number of components-an alternative to reversible jump methods[END_REF] proposed a continuous time version of RJMCMC, based on earlier ideas of [START_REF] Geyer | Simulation procedures and likelihood inference for point processes[END_REF], but with similar properties [START_REF] Cappé | Reversible jump, birth-and-death, and more general continuous time MCMC samplers[END_REF], while [START_REF] Brooks | Efficient construction of reversible jump Markov chain Monte Carlo proposal distributions (with discussion)[END_REF] made proposals for increasing the efficiency of the moves.

In retrospect, while reversible jump is somehow unavoidable in the processing of very large numbers of models under comparison, as for instance in variable selection [START_REF] Marin | Bayesian Core[END_REF], the implementation of a complex algorithm like RJMCMC for the comparison of a few models is somewhat of an overkill since there exist alternative solutions based on model specific MCMC chains, for example [START_REF] Chen | Monte Carlo Methods in Bayesian Computation[END_REF].

Regeneration and the CLT

The Ergodic Theorem (see, for example, Robert and Casella 2004, Theorem 6.63) is essentially the Strong Law of Large Numbers rewritten for Markov chains. If X 1 , X 2 , • • • X n is a Markov chain with stationary distribution π, and h(•) is a function with finite variance, then under fairly mild conditions,

lim n→∞ hn = h(x)π(x) dx = E π h(X), (3)
almost everywhere, where hn = (1/n) n i=1 h(X i ) . To monitor this convergence, we would like to appeal to a Central Limit Theorem (CLT) and use the fact that √ n( hn -E π h(X))

Varh(X) → N(0, 1), ( 4) but there are two roadblocks to this. First, convergence to normality is strongly affected by the lack of independence. To get CLTs for Markov chains, we can use a result of [START_REF] Kipnis | Central limit theorem for additive functionals of reversible Markov processes and applications to simple exclusions[END_REF], which requires the chain to be reversible (a fact that holds for Metropolis-Hastings chains), or we must delve into "mixing conditions" (Billingsley 1995, Section 27), which are typically not easy to verify. However, [START_REF] Chan | Discussion of "Markov chains for exploring posterior distribution[END_REF] showed how the condition of geometric ergodicity could be used to establish CLTs for Markov chains. But getting the convergence is only half of the problem. In order to use (4), we must be able to consistently estimate the variance, which turns out to be another difficult endeavor. The "naïve" estimate of the usual standard error is not consistent in the dependent case (try the simple calculation where the X i are equicorrelated), and the most promising paths for consistent variance estimates seems to be through regeneration and batch means.

The theory of regeneration uses the concept of a split chain (Athreya andNey 1978, Robert andCasella 2004, Chapter 6), and allows us to independently restart the chain while preserving the stationary distribution. These independent "tours" then allow the calculation of consistent variance estimates and honest monitoring of convergence through (4). Early work on applying regeneration to MCMC chains was done by [START_REF] Mykland | Regeneration in Markov chain samplers[END_REF] and [START_REF] Robert | Convergence control techniques for MCMC algorithms[END_REF], who showed how to construct the chains and use them for variance calculations and diagnostics (see also [START_REF] Guihenneuc-Jouyaux | Finite Markov chain convergence results and MCMC convergence assessments[END_REF], as well as deriving adaptive MCMC algorithms [START_REF] Gilks | Adaptive Markov chain Monte Carlo[END_REF]. [START_REF] Rosenthal | Minorization conditions and convergence rates for Markov chain Monte Carlo[END_REF] also showed how to construct and use regenerative chains, and much of this work is reviewed in [START_REF] Jones | Honest exploration of intractable probability distributions via Markov chain Monte Carlo[END_REF]. The most interesting and practical developments, however, are in [START_REF] Hobert | On the applicability of regenerative simulation in Markov chain Monte Carlo[END_REF] and [START_REF] Jones | Fixed-width output analysis for markov chain monte carlo[END_REF], where consistent estimators are constructed for Varh(X), allowing valid monitoring of convergence in chains that satisfy the CLT. Interestingly, although [START_REF] Hobert | On the applicability of regenerative simulation in Markov chain Monte Carlo[END_REF] uses regeneration, [START_REF] Jones | Fixed-width output analysis for markov chain monte carlo[END_REF] get their consistent estimators thorough another technique, that of cumulative batch means.

Conclusion

The impact of Gibbs sampling and MCMC was to, almost instantaneously, change our entire method of thinking and attacking problems, representing a paradigm shift in the words of the historian of science Thomas Kuhn [START_REF] Kuhn | The Structure of scientific Revolutions, Third Edition[END_REF]. Now, the collection of real problems that we could solve grew almost without bound. Markov chain Monte Carlo changed our emphasis from "closed form" solutions to algorithms, expanded our impact to solving "real" applied problems, expanded our impact to improving numerical algorithms using statistical ideas, and led us into a world where "exact" now means "simulated"! This has truly been a quantum leap in the evolution of the field of statistics, and the evidence is that there are no signs of slowing down. Although the "explosion" is over, the current work is going deeper into theory and applications, and continues to expand our horizons and influence by increasing our ability to solve even bigger and more important problems. The size of the data sets, and of the models (for example in genomics or climatology) is something that could not have been conceived 60 years ago, when Ulam and von

Neumann invented the Monte Carlo method. Now we continue to plod on, and hope that the advances that we make here will, in some way, help our colleagues 60 years in the future solve the problems that we cannot yet conceive!

A Appendix: Workshop on Bayesian Computation

  , z t,k ) , z t,k ∼ πt-1 (z|x) , k = 1, . . . , m, that is, by replicating m times the simulations from the current approximation πt-1 (z|x) of the marginal posterior distribution of the missing data. This focus on estimation of the posterior distribution connected the original Data Augmentation algorithm to EM, as pointed out by Dempster in the discussion. Although the discussion by Morris gets very close to the two-stage Gibbs sampler for hierarchical models, he is still concerned about doing m iterations, and worries about how costly that would be. Tanner and Wong mention taking m = 1 at the end of the paper, referring to this as an "extreme case".In a sense,[START_REF] Tanner | The calculation of posterior distributions by data augmentation[END_REF] was still too close to Rubin's 1978 multiple imputation to start a (new) revolution. Yet another reason for this may be that the theoretical background was based on functional analysis rather than Markov chain theory, which needed, in particular, for the Markov kernel to be uniformly bounded and equicontinuous. This may have discouraged potential applicants as requiring too much math! The authors of this review were fortunate enough to attend many focused conferences during this time, where we were able to witness the explosion of Gibbs sampling. In the summer of 1986 in Bowling Green, Ohio, Adrian Smith gave a series of ten lectures on hierarchical models. Although there was a lot of computing mentioned, the Gibbs sampler was not fully developed yet. In another lecture by Adrian Smith in June 1989 at a Bayesian workshop in Sherbrooke, Québec, Adrian revealed for the first time (?) the generic features of Gibbs sampler, and we still remember vividly the shock induced on ourselves and on the whole audience by the sheer breadth of the method! This development of Gibbs sampling, MCMC, and the resulting seminal paper of Gelfand andSmith (1990) was an epiphany in the world of Statistics.Definition: epiphany n. A spiritual event in which the essence of a given object of manifestation appears to the subject, as in a sudden flash of recognition.The explosion had begun, and just two years later, at an MCMC conference at Ohio State University organized by Alan Gelfand, Prem Goel, and Adrian Smith, there were three full days of talks. The presenters at the conference read like a Who's Who of MCMC, and the level, intensity and impact of that conference, and the subsequent research, is immeasurable.The program of the conference is reproduced in Appendix A. Approximately one year later, in May of 1992, there was a meeting of the Royal Statistical Society on "The Gibbs sampler and other Markov chain Monte Carlo methods", where four papers were presented followed by much discussion. The papers appear in the first volume of JRSSB in 1993, together with 49 (!) pages of discussion, again by the Who's Who of MCMC, and the excitement is clearly evident in the writings (even though the theory and implementation were not always perfectly understood).

  a) Theoretical Aspect of Iterative Sampling, Chair: Adrian Smith 1) Martin Tanner, University of Rochester: EM, MCEM, DA and PMDA 2) Nick Polson, Carnegie Mellon University: On the Convergence of the Gibbs Sampler and its Rate 3) Wing-Hung Wong, Augustin Kong, and Jun Liu, University of Chicago: Correlation Structure and Convergence of the Gibbs Sampler and Related Algorithms (b) Applications -I, Chair: Prem Goel 1) Nick Lange, Brown University, Brad Carlin, Carnegie Mellon University and Alan Gelfand, University of Connecticut : Hierarchical Bayes Models for Progression of HIV Infection 2) Cliff Litton, Nottingham University, England: Archaeological Applications of Gibbs Sampling 3) Jonas Mockus, Lithuanian Academy of Sciences, Vilnius: Bayesian Approach to Global and Stochastic Optimization • Saturday, Feb. 16, 1991 (a) Posterior Simulation and Markov Sampling, Chair: Alan Gelfand 1) Luke Tierney, University of Minnesota: Exploring Posterior Distributions Using Markov Chains 2) Peter Mueller, Purdue University: A Generic Approach to Posterior Integration and Bayesian Sampling 3) Andrew Gelman, University of California, Berkeley and Donald P. Rubin, Harvard University: On the Routine Use of Markov Chains for Simulations 4) Jon Wakefield, Imperial College, London: Parameterization Issues in Gibbs Sampling 5) Panickos Palettas, Virginia Polytechnic Institute: Acceptance-Rejection Method in Posterior Computations (b) Applications -II, Chair: Mark Berliner 1) David Stephens, Imperial College, London: Gene Mapping Via Gibbs Sampling 2) Constantine Gatsonis, Harvard University: Random Effects Model for Ordinal Categorical Data with an application to ROC Analysis 3) Arnold Zellner, University of Chicago, Luc Bauwens, Université de Louvainla-Neuve, and Herman Van Dijk, Rotterdam University: Bayesian Specification Analysis and Estimation of Simultaneous Equation Models using Monte Carlo Methods (c) Adaptive Sampling, Chair: Carl Morris 1) Mike Evans, University of Toronto and Carnegie Mellon University: Some Uses of Adaptive Importance Sampling and Chaining 2) Wally Gilks, Medical Research Council, Cambridge, England: Adaptive Rejection Sampling 3) Mike West, Duke University: Mixture Model Approximations, Sequential Updating and Dynamic Models • Sunday, Feb. 17, 1991 (a) Generalized Linear and Nonlinear Models, Chair: Rob Kass 1) Ruey Tsay, and Robert McCulloch, University of Chicago: Bayesian Analysis of Autoregressive Time Series 2) Christian Ritter, University of Wisconsin: Sampling Based Inference in Non Linear Regression 3) William DuMouchel, BBN Software, Boston: Application of the Gibbs Sampler to Variance Component Modeling 4) James Albert, Bowling Green University and Sidhartha Chib, Washington University, St. Louis: Bayesian Regression Analysis of Binary Data 5) Edwin Green and William Strawderman, Rutgers University: Bayes Estimates for the Linear Model with Unequal Variances (b) Maximum Likelihood and Weighted Bootstrapping, Chair: George Casella 1) Adrian Raftery, and Michael Newton, University of Washington: Approximate Bayesian Inference by the Weighted Bootstrap 2) Charles Geyer, Universlty of Chicago: Monte Carlo Maximum Likelihood via Gibbs Sampling 3) Elizabeth Thompson, University of Washington: Stochastic Simulation for Complex Genetic Analysis (c) Panel Discussion -Future of Bayesian Inference using Stochastic Simulation, Chair: Prem Gael • Panel -Jim Berger, Alan Gelfand, and Adrian Smith

  This section contains the program of the Workshop on Bayesian Computation via Stochastic Simulation, held at Ohio State University, February 15-17, 1991. The organizers were Alan Gelfand, University of Connecticut, Prem Goel, Ohio State University, and Adrian Smith,

	Imperial College, London.

•

Friday, Feb. 15, 1991

The atomic bomb construction did not involve simulation techniques, even though the subsequent development of the hydrogen bomb did

.3 MANIAC stands for Mathematical Analyzer, Numerical Integrator and Computer.

On a somber note, Edward Teller later testified against Robert Oppenheimer in the McCarthy trials and, much later, was a fervent proponent of the "Star Wars" defense system under the Reagan administration.

In fact, Hastings starts by mentioning a decomposition of the target distribution into a product of one-dimensional conditional distributions but this falls short of an early Gibbs sampler!

[START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their application[END_REF] is one of the ten papers reproduced in the Biometrika 100th anniversary volume by[START_REF] Titterington | Biometrika: One Hundred Years[END_REF].

project ANR=05-BLAN-0299 Adap'MC.

ported by National Science Foundation Grants DMS-04-05543, DMS-0631632 and SES-0631588.


Maryam Khalili

Thierry Candresse

Chantal Faure

Armelle Marais
email: armelle.marais-colombel@inrae.fr

Complete genome sequence of almond luteovirus 1, a novel luteovirus species infecting almond

Keywords: stone fruits, virome, high-throughput sequencing, Luteoviridae, Prunus Arch. Virol

In this study, we report the complete genome sequence of a novel luteovirus detected in almonds using high-throughput sequencing. The genome of the new luteovirus comprises 5,047 nucleotides, and its genomic organization is similar to that of the recently described nectarine stem pitting associated virus (NSPaV), with only four open reading frames encoding for replication related proteins, the coat protein (CP) and a CP readthrough protein involved in the aphid transmission of luteovirids. Phylogenic and pairwise distance analyses showed that this virus shares 79% and 57.8% amino acid identities in the P1-P2 fusion protein and the P3-P5 protein, respectively, with the most closely related luteovirus NSPaV, suggesting that it represents a novel species for which the name of Almond associated luteovirus 1 is proposed.

To our knowledge, this is the first report of an almond-infecting luteovirus.

INTRODUCTION

Almond (Prunus dulcis) is believed to have originated from Southwest Asia where it was initially brought into domestication and is today cultivated in many countries. Grafttransmissible pathogens (viruses, viroids and phytoplasmas) cause economically important diseases in almond, in particular viruses of the genera Ilarvirus, Trichovirus, Nepovirus and Potyvirus [START_REF] Savino | Virus and virus-like diseases of pome and stone fruits[END_REF]. However, in contrast to several other Prunus species, no luteovirus has thus far been reported from almond.

Luteoviruses belong to the family Luteoviridae, which also comprises the genera Enamovirus and Polerovirus [START_REF] Domier | Family Luteovirideae[END_REF]. Three luteoviruses, nectarine stem pitting-associated virus (NSPaV) [START_REF] Bag | Detection of a New Luteovirus in Imported Nectarine Trees: A Case Study to Propose Adoption of Metagenomics in Post-Entry Quarantine[END_REF], cherry-associated luteovirus (ChALV) [START_REF] Lenz | Identification and characterization of a new member of the genus Luteovirus from cherry[END_REF] and peach-associated luteovirus (PaLV) [START_REF] Wu | Molecular characterization of a novel luteovirus from peach identified by high-throughput sequencing[END_REF] have recently been discovered in Prunus. The genome of luteoviruses varies between 5 and 5.9 kilobases (kb) [START_REF] Domier | Family Luteovirideae[END_REF]. Among luteoviruses described so far, NSPaV is the sole member with only four open reading frames (ORFs) [START_REF] Bag | Detection of a New Luteovirus in Imported Nectarine Trees: A Case Study to Propose Adoption of Metagenomics in Post-Entry Quarantine[END_REF], whereas other members contain six or seven ORFs [START_REF] Domier | Family Luteovirideae[END_REF].

As part of a virus discovery effort, double-stranded RNAs were purified from the leaves of an almond tree without obvious viral symptoms collected in Turkey in the spring of 2019.

Following reverse transcription and random amplification [START_REF] Marais | Viral Double-Stranded RNAs (dsRNAs) from Plants: Alternative Nucleic Acid Substrates for High-Throughput Sequencing[END_REF], the obtained cDNAs were sequenced (Hiseq3000 Illumina platform, 2x150 bp). After reads cleaning, assembly into contigs and Blast-based annotation [START_REF] Candresse | Association of Little cherry virus 1 (LChV1) with the Shirofungen stunt disease and characterization of the genome of a divergent LChV1 isolate[END_REF], a single 5 kb contig showing significant nucleotide (nt) identity with luteoviruses was identified. This contig integrated 35,411 reads (34.3% of total reads) for an average coverage of 782x and showed 71% nt identity with the NSPaV genome (NC027211) which spanned from nt 39 to 4,991. The genome sequence was completed by determining the genome ends using the SMARTer RACE Kit (Takara) and has been deposited in GenBank under accession number MT362517.

The complete genome of the almond luteovirus is 5,047 nt long and comprises four ORFs, homologous to the corresponding NSPaV ORFs (Fig. 1). ORF1 encodes the putative RNAdependent RNA polymerase (RdRp) of 329 amino acids (aa), while an 848 aa fusion protein is expressed by a -1 ribosomal frameshift between ORF1 and ORF2. ORF3 encodes the putative 204 aa coat protein (CP) and ORF5 is expressed as a readthrough fusion to the CP (550 aa) following the suppression of the CP stop codon (Fig. 1). With 637 nt, the long 3' untranslated region (UTR) has a length similar to that of NSPaV, while the 5' UTR is 149 nt long. Similar to NSPaV, ORF5 (345 aa) was found to be significantly shorter than in other luteoviruses, and neither ORF3a [START_REF] Smirnova | Discovery of a Small Non-AUG-Initiated ORF in Poleroviruses and Luteoviruses That Is Required for Long-Distance Movement[END_REF] nor ORF4 (putative movement protein) nor ORF6 could be identified by scanning the genome and analyzing translation in all reading frames. A phylogenetic analysis using complete luteovirus genomes (GenBank) confirms that the almond virus belongs to the genus Luteovirus and has the closest affinities with NSPaV with which it shares 71.5% nt identity over the whole genome (Supplementary Fig. 1). Phylogenetic trees reconstructed using amino acid sequences of the RdRp or CP showed that NSPaV and the almond virus form a specific branch, the RdRp tree showing a lower branching as compared to the CP tree (Fig. 2A andB). Pairwise comparisons for the various viral proteins showed that the almond virus P1-P2 fusion protein shares 49.6%-79% aa identity with those of luteoviruses, the highest being NSPaV. For the CP-P5 readthrough protein, the level of identity was even lower (27.2%-57.8%). Given the currently accepted molecular species demarcation criteria for luteoviruses (more than 10% in aa divergence in any gene product [START_REF] Domier | Family Luteovirideae[END_REF]), the almond virus defines a novel Luteovirus species for which the name Almond luteovirus 1 (AlLV1) is proposed. In order to clarify potential recombination events with related viruses, a recombination analysis was performed using RDP 4 on a multiple alignment of representative isolates of recognized luteovirus species but failed to identify statistically supported recombination events involving AlLV1.

To estimate the incidence of AILV1 in Prunus, detection primers were designed (NLuteo-PolF1 5'-GCACATGTTTCGAGGAGATG-3', NLuteo-PolR1 5'-CCGTCCTCAGATCTTTCCAT -3') targeting a fragment of 505 nt in the polymerase gene. Total nucleic acids were extracted [START_REF] Foissac | Polyvalent degenerate oligonucleotides reverse transcription-polymerase chain reaction: a polyvalent detection and characterization tool for trichoviruses, capilloviruses, and foveaviruses[END_REF] from 105 Prunus samples (P. dulcis, P. sibirica, P. avium, P. armeniaca, P. cerasifera, P. domestica, P. mume and P. persica) collected from Azerbaijan, Kazakhstan, China, Italy and France and submitted to a two-step RT-PCR [START_REF] Marais | Development of a polyvalent RT-PCR detection assay covering the genetic diversity of Cherry capillovirus A[END_REF]. Besides the positive sample, none of the tested samples were found to be infected by AlLV1. To our best knowledge, AILV1 is the first described luteovirus infecting almonds. Further studies are required to estimate its distribution in stone fruits, as well as its potential pathogenicity. As for NSPaV, the absence of ORFs 3a, 4 and 6 and the presence of a shorter readthrough domain might suggest that its biology could significantly differ from that of other luteoviruses.

Figure1.

  Figure1. Schematic representation of the genome organization of almond luteovirus 1.Boxes represent the Open Reading Frames (ORF) and numbers indicate the coordinates of

Figure 2 .

 2 Figure 2. Neighbor-joining phylogenetic trees reconstructed using the alignments of

Acknowledgements

The authors thank the INRAE GenoToul Platform (Toulouse, France) for HTS, and Eden Darnige for proofreading the manuscript.

Compliance with ethical standards: this manuscript presents original work that does not involve studies with humans or animals.

Conflict of interest: the authors declare that they have no conflict of interest.

Supplementary


Danesh Sadjad

Mesgaran

Juan Haladjian

Stefan Nüske

Dorothée Ledoux

David Humphries

Lene Munksgaard

Isabelle Veissier

Lameness detection and scoring

Keywords:

Introduction

Enting [START_REF] Enting | Economic losses due to clinical lameness in dairy cattle[END_REF] concluded lameness, from an economic perspective, as the third most costly health disease, following mastitis and reproductive failure issues, in cattle units. Archer [START_REF] Archer | Lameness in UK dairy cows: a review of the current status[END_REF] estimated the incidence rate of lameness in the United Kingdom cattle herds roughly 50 cases/100 cows in a year; nevertheless, due to poor correlation between incidence rates and records of treatments in farms [START_REF] Whay | Assessment of the welfare of dairy cattle using animal-based measurements: direct observations and investigation of farm records[END_REF], the actual number seems to be higher. Surprisingly, the significance of lameness associated with cattle welfare, health and profitability of the unit has been greatly underestimated [START_REF] Van Nuffel | Lameness detection in dairy cows: part 1. How to distinguish between non-lame and lame cows based on differences in locomotion or behavior[END_REF]. Recent works have shown a clear link between BCS and hook condition of cows with the development of lameness in these animals [START_REF] Randall | Low body condition predisposes cattle to lameness: an 8-year study of one dairy herd[END_REF]. Lameness is a multifactorial and progressive issue where different detriments contribute to its development via complex interactions. Detection of lame cattle can be facilitated through description of the animals' gait characteristics in a numerical scaling system known as locomotion scoring [START_REF] Sprecher | A lameness scoring system that uses posture and gait to predict dairy cattle reproductive performance[END_REF]. The total number of visual (manual) locomotion scoring systems can reach up to 25, where differences lie mostly in the used scales, characterization of cows' gait, and posture [START_REF] Schlageter-Tello | Effect of merging levels of locomotion scores for dairy cows on intra-and interrater reliability and agreement[END_REF]. Automated locomotion scoring tools would be a big advantage for regular monitoring of lameness in the herd. Three methods that are commonly engaged with automated systems are: kinetic, kinematic and indirect. The kinetic and kinematic approaches measure the forces, involved in locomotion, and time and distance of variables, associated to limb movement, respectively [START_REF] Flower | Gait assessment in dairy cattle[END_REF]. The indirect method simply exploits behavioral or production data as indicators for impaired locomotion. The automated tools/instruments that will be developed, based upon either of the aforementioned approaches, need to be validated with a 'reference' method. This usually is done by comparing with manual scoring; however, it is noteworthy that manual scoring systems have their own set of limitations [START_REF] Chapinal | Using gait score, walking speed, and lying behavior to detect hoof lesions in dairy cows[END_REF].

Prerequisites

This guideline details steps for detecting and assessing lameness in cattle within the experimental units. These steps are written, presuming that lameness detection is carried out in the units either manually (visual detection) or via automatic lameness monitoring systems. The Animal Trait Ontology for Livestock (ATOL) numbers linked with this guideline are: ATOL_0000835 and ATOL_0001780 (for complete list of ATOL, please visit https://www.atol-ontology.com/en/erter-2/).

A -Lameness detection based on visual evaluation

1. The established scoring system within the experimental units should be easy to apply routinely and be time efficient. As recommending a distinct scoring system which is suitable for all experimental units is difficult, use should be limited to just a few systems such as the 3-point scale, e.g. Welfare Quality protocol (http://www.welfarequality.net/en-us/reports) or the 4-point scaling system as recommended by AHDB (http://www.welfarequality.net/media/1121/wqr11.pdf).

Regardless of the type of scoring system, it is recommended that the methods used for visual examination are justified according to the literature or in-field practice. 2. A systematic lameness assessment of all cows should be carried out every 3-4 weeks (unless the experiment dictates another frequency), whilst any new case of lameness detected in between is also recorded. 3. Scoring of dairy cattle should be done after milking, as animals return to the field or housing. 4. Operators in the unit should be aware of the close association between lameness and longer lying times and standing in passageway, as well as altered feeding (grazing) behaviour (unless the experiment design assesses other factors). Dairy cows in the unit that display signs of lameness require careful monitoring. 5. The surface, animals walk over during lameness scoring, should have a flat and even surface, and should not be slippery or muddy as this will impact gait and lead to inaccurate locomotion scoring. 6. Animals should always be evaluated in the same location to reduce variation in scores for repeat observations. 7. Ideally, the same operator should locomotion-score the animals from a given herd each time to reduce observer error. If there are several operators, their consistency must be checked regularly (twice a year is recommended). Regular calibration between observers should be carried out. 8. The operator should not induce any stress to the animal during the scoring procedure. They should not make any sudden movements or force animals to move forwards. Scoring agitated cows will produce inaccurate scores. 9. The operator responsible for locomotion scoring should be trained and retrained by an experienced observer, familiar with the scoring system utilised. This will significantly increase the agreement level between and within observers.

B -Automated lameness detection systems

1. Data from automated systems based on visual characteristics of the animal, e.g. change in locomotion or behaviour, must be validated with manual locomotion scoring as a reference ('gold standard'). 2. The automated system should be calibrated according to the reference visual method as stated in point 1 of section A . Users are highly recommended to see the guideline for sensor validation within this chapter of the book. 3. Lameness detection with automated tools should be conducted during the normal daily routine (i.e. after leaving the milking parlour) of cows/heifers and ideally without involvement of the operators/staff members (e.g., the platform with loading cells should be placed behind the milking parlor or at the exit to the pasture during grazing season). Animals forced towards an experimental set-up, may hide their weakness or pain due to their stoic nature [START_REF] Weary | Board-invited review: Using behavior to predict and identify ill health in animals[END_REF]. 4. The operator should ensure that dairy cows/heifers do not block the area around the measurement zone when walk-over devices, such as StepMetrix, are being used as this will produce faulty data.

C -Lameness detection using a wearable sensor

The present guideline details steps for collection and data analysis for lameness detection applications. This document has been written, based on a recently developed automated lameness detection system using a wearable motion sensor [START_REF] Haladjian | A Wearable Sensor System for Lameness Detection in Dairy Cattle[END_REF]. The developed system made use of the MicroHub (Interactive wear, Starnberg, Germany). Where the research unit intends to utilise the MicroHub for data collection, a mobile phone is necessary to start and stop data recordings. A detailed description of this lameness detection system has been created by Haladjian [START_REF] Haladjian | A Wearable Sensor System for Lameness Detection in Dairy Cattle[END_REF]. The Animal Trait Ontology for Livestock (ATOL) number linked with this guideline is: ATOL_0000351, AHOL_0005001, EOl_0000001, EOL_0001606, EOL_0001706 and EOL_0001712 (for complete list of ATOL, please visit https://www.atolontology.com/en/erter-2/).

Preparing for data collection:

The data collection should be planned in advance, based on the particular purpose of the study, and documented in a collection plan. A data collection plan should include: a) list of materials to be used; b) concise instructions to be followed by the person in charge of the data collection (e.g., how to attach the sensor device, how to orient it, how to start recording data); c) list of actions to be followed during the experiment (e.g., list of cows to be recorded, what needs to be recorded ,and for how long). 2. The selection of the cows that will participate in the experiment should be done together with the veterinarian responsible for the cows. 3. Modern machine learning algorithms achieve better performances when trained with larger amounts of data. However, cows should not be forced to walk for longer than 10 min per session not only for ethical reasons, but also to avoid a change in the walking pattern due to fatigue. It is recommended that data collection sessions are spread over more than one day for a cow. 4. The batteries of devices (e.g., mobile phone, sensor devices) used during the experiment should be fully charged. 5. It is recommended that the operator conducts a short (e.g. 10s-long) data collection to ensure proper function of the data collection system (e.g. mobile device and sensor device) and to reduce the chance of human errors before attaching the sensor to a cow. 6. Unless the experimental design dictates otherwise, the newly established detection system should be tested on a wide range of animals, based on their physiological conditions (parity, weight, breed, and lactation stage), environmental condition (indoor housing, grazing management), or external parameters (time of the year or day). 7. The health condition of the cows selected for data collection should be considered. Cows that are pregnant, undergoing oestrus, or obviously sick (e.g. lame) should be excluded from the experiment. 8. Data collection: Regardless of the type of motion sensor used in the unit, the operator should not touch the data collecting device or recording instruments with wet hands. Data collection can be conducted by a single operator if they have experience manipulating cows. However, it is useful to have a second individual to help move cows and restrain them to attach or remove the device. 9. The animal must be restrained so that the sensor device can be firmly attached. The optimal place for attachment of the device must be determined, based on the particular application. For lameness, it is recommended to attach the device to the hind left limb as lameness affects mostly hind limbs. How the device is oriented and to which hind limb it is attached does usually not matter, as long as these parameters are consistent across cows and recordings. Figure 1 shows a device attached to a cow and the recommended orientation of the device. . The animal should be encouraged to walk to allow analysis of its gait. The goal of the data collection is to collect data that resembles the context in which the lameness detection system would be used. When encouraging a cow to walk, the operator should therefore maintain enough distance to the cow, so that it can walk at its own pace. Causing fear or stress to the animal could easily bias the algorithm.

11. Researchers should pay attention to the cow's gaze, particularly when crossing the milking parlour and free pens. Cows will try to enter the milking parlour or free pens to avoid having to participate in a data collection experiment. If possible, the environment should be arranged in such a way that cows walk in a part of the stable that has no access to the milking parlour. 12. Operators should document when a cow stops walking for a long time, so that this event can be considered when studying the data afterwards. If a cow enters a pen, it can be encouraged to leave it simply by waving. It is desirable to generate a data set with a similar amount of data across cows. However, only data recorded while a cow walks is useful for lameness detection. Therefore, the data collection should last longer for cows that entered the milking robot or stopped walking for a long period. 13. It is often necessary to collect reference data to assess and annotate collected motion data. Some studies may use manual annotation, noting whether a cow walked normally or with a plastic block attached to a hind hoof as a reference. Other studies may require reference values for the entire data set. For example, if a lameness detection system counts the number of strides performed by a cow, reference values for the data would be needed to indicate at which moment in time a cow performed a stride. A common reference data source is video. When using reference data, the motion and reference data sets should be synchronised. This can be done by identifying and matching specific events in the recorded motion and reference data. 14. MicroHub: The MicroHub should be turned on by pressing and holding the button on it for 5 s.

The light will blink, then display a steady green LED when the MicroHub is turned on. On the app, tapping on the 'scan' button will scan for devices nearby with Bluetooth Low Energy (BLE) capability (other BLE devices may show in addition of the MicroHub). From the devices shown, MicroHub should be selected and connected to by tapping 'connect'. When connected, the MicroHub displays a blue LED. 15. Use the 'start/stop' button on the app to start recording. When the MicroHub start recording data, it displays a red LED and the app displays the text: 'Logging Since: sec' Once recording, tapping again on the 'start/stop' button on the app will cause the MicroHub to stop recording data. Every time, a recording is started, a new file is created in the MicroHub and stored on its SD card. 16. Once recording has started, disconnecting the app from the MicroHub will save its battery. When disconnected, the blue LED should turn off. The MicroHub can still record data when it is not connected to the app. 17. Before the SD card can be safely removed from the MicroHub, the MicroHub should be turned off. It can be turned off by pressing the white button for up to 4 s. The SD card can be taken out by pushing, then releasing it, and then removing. Data on the miniSD card can be imported to a computer with an adapter. The SD card should contain a text file, named according to the date and time when the recording started. For example, if the recording was started on April 26th 2018 at 14:03 h, the file will be named '2018.04.26_14.03.50.969.txt'. The file name should be assessed against the data collection plan and notes taken during the data collection. 18. Data analysis: Once data has been collected, an initial check should be carried out to ensure the recording was successful. It is advisable to import the data immediately following data collection, whilst researchers remember the flow of the data collection session (e.g. the order of cows if anything differed from the recording plan). A mistake at this stage could invalidate the data analysis. 19. Several things need to be considered when preparing data for analysis: The motion tracker may store data in text format. This is convenient for data collection purposes (e.g. to quickly assess whether data stored in a file is properly stored). However, loading the entire data set in a text format can slow down the data analysis process, as files must be loaded several times. Therefore, data should be stored in a binary format, which will occupy considerably less memory and therefore load faster. 20. If Matlab is used for the data analysis, variables in memory can be saved with a simple 'save<' command. Before starting data analysis, data should be cleaned up by removing unwanted sections. For example, if the sensor is attached to the cow after starting the recording, data corresponding to the action of attaching it cow should be removed from the analysis. 21. Specific gestures (e.g. shaking the sensor), used to synchronise the sensor data to a reference data set, should also be removed. 22. Rename data files to uniquely identify a cow and the kind of data it contains in a structured way, so that the different files can be loaded automatically in the script used for analysis. For example 'cowX_Y' where X indicates the ID number of the cow and Y the type of recording done for the particular cow (e.g. left to indicate the cow walked with a block attached to its left hoof). 23. Researchers should manually test and visualise data processing methods on a small fraction of data firstly. However, to obtain a consistent and statistically significant metric to compare different data processing methods, a metric should be computed automatically for the entire data set. This metric might be the classification accuracy or the F1-score. 24. System design: These steps are meant for long-term deployment of a detection system, regardless of the type of sensor utilised. Design decisions should be made, regarding the sensor device, computations that will be executed on the device, and the architecture of the system. 25. Energy consumption of the device should be minimised to avoid frequent battery replacements.

For this, computation should be done locally on the wearable device and the computed results transmitted, at as sparse time intervals as acceptable for the particular application. 26. The sensor should be designed to minimise the risk that the animal injures itself, other cows in the stable (e.g. by knocking the device into them), or the device itself. 27. The system should be designed to precisely generate an individual gait profile for each cow and detect changes in their locomotion. 28. The system's data receiver should be installed at the milking parlour. Cattle usually visit the milking parlour at least twice a day and remain still in the same area for several minutes during milking. This period represents an ideal opportunity for the data transmission.

Figure 1 :

 1 Figure 1: Sensor device attached to a left hind leg and recommended axis-orientation of the wearable device

Methods in cattle physiology and behaviour research -Recommendations from the SmartCow consortium

Copyright: © 2020 Sadjad Danesh Mesgaran et al. This is an Open Access publication distributed under the terms of the Creative Commons Attribution 4.0 International License. See license information at https://creativecommons.org/licenses/by/4.0/


Quieter and more civilised but not more accountable

In the face of the Covid-19 crisis, the Hansard Society wonders, in one of its last April publications, "how parliamentarians can conduct their core constitutional duties of holding the government to account, when [they] are all required to adopt rigorous social distancing".

This question is even more relevant in relation to PMQs insofar as it is the most attended and the most watched parliamentary event of the week.

Due to social distancing far fewer MPs have been able to attend PMQs, and that resulted in less noise being made and fewer raucous attacks levelled (by the parties).

Sky News John Craig recommends that "[they] need the Commons to return to its usual boisterous, raucous and spontaneous normality if government ministers are to be held truly accountable".

Those remarks beg the question of how disorder is an effective instrument in holding the government to account. And what "truly accountable" means in the specific context of PMQs.

The interruption as a tool to hold the government to account I'd like to focus on two types of disorderly behaviour : interruptions on the one hand, and on the other hand : unparliamentary language -which is here defined as : threats and attacks to the face 1 (FTAs like insults, soundbites, offending remarks…).

Besides the adversarial and confrontational style of parliamentary politics, what may account for the extensive use of FTAs (disorderly language) lies at the core of the function of PMQs : holding the gvt to account.

According to Graham P. Thomas's 2 definition of accountability, "criticising the government's failure, incompetence or deceit" are FTAs per se, which inevitably causes disturbance especially when the PM counter-attacks to save face. More often than not, those departures from the rules are tolerated by the Speaker, very few FTAs are sanctionned ; Sandra Harris 3 mentions a very high level of tolerance of 'illegal interventions'. This tolerance towards disorderly language lies at the root of accountability. So it is no wonder that MPs who want to make their point stronger should strike harder, which inevitably leads to an increase in the number and in the intensity of the attacks. In other words, the level of expected accountability might be reflected in the intensity and frequency of the attacks.

Whether they are challenging or dismissive backchanneling interruptions (if we use Cornilia Ilie's typology 4 ), they serve the purpose of interrupting the scripted speech of the PM urging him/her to answer more spontaneously Those impromptu answers are valuable in apparaising the government's accountability. The PM will be judged more on the answers he will provide off the cuff than on his ready-made answers. The consequence (or the purpose) of such interruptions is to warn people about a message that is being delivered ; they function as attention getters. Whether the interruption points out an actual inaccuracy or even a downright lie, it doesn't matter. The government's accountability doesn't relate so much to the message or to the speech of the PM but to those signals sent out to the pubic prompting them to call into question specific points or even to a greater extent call into question the person's credibility.

What does "true accountability" mean in the context of PMQs ?

MPs represent their constituents in Parliament, and defend their interests. Unlike other parliamentary sessions, PMQs is a vehicle for individual backbenchers to raise individual grievances of constituents. MPs are the voice and in a way the body of their constituents in Parliament. In this theatrical performance5 backbenchers do not remain a passive audience. They make themselves and the people's voice heard by causing disturbance, by interrupting the speaking debaters. Backbenchers who are gathered in the House acting on behalf of the citizens take part in the debate through the noise they make in reaction to the speech of the leaders. This way of interacting with the executive in order to hold them accountable is popular in essence. By popular I mean that it is the most watched parliamentary event of the week and that also means that PMQs involves a way of discurive interaction and styles of communication pertaining to and deriving from the people.

In ancient Athens, the thorubos (meaning the noise of the participants) was the expression of the sovereignty of the people (dèmos). Making noise and interrupting the debaters represented an active participation in the debate. The disruptions and roar of the backbenchers is the popular expression of the distrust, the discontent or the indignation of the people and in that case those feelings are not mediated by long-winded technical pre-scripted speeches.

Unquestionably, this is a more direct way of hoding the government to account.

Even though the House debates have been subject to a new configuration where the paradigms in terms of disorderly behaviour have critically changed, MPs have been reluctant (unable ?) to change their discursive model and are still prone to disorderly attacks and interruptions. It took MPs only seven sessions with social distancing before going back to their old ways.

It seems that causing some kind of disorder through FTAs and interruptions is a natural and inescapable process of holding the governmentt to account, and that is what 'true accountability' implies, that is to say making the voice and the noise of the people heard by using the mode of communication of the ordinary citizens.

That seems to suggest that the distance between the government and the people should be reduced to increase government's accountability.

How is the distance reduced ?

In the orderly question-answer sequences those disruptions act as metadiscursive (or rather metadramatic I should say) short-cuts or substitutes for longer, more elaborate answers or counterarguments. This simplification of the message makes it easier for people to grasp the meaning of the spontaneous intervention of a member.

For Cornelia Ilie 6 the emotional force of insults that is correlated to the (over)-simplification of the message is tremendous and outweighs their rational force. I would go as far as contending that insults produce the same effetcs as symbols. They condense a wide range of meanings and emotions.

For Murray Edelman 7 , people think in terms of stereotypes and oversimplifications due to some incapacity to "tolerate ambiguous and complex situations and respond chiefly to symbols that oversimplify and distort". There is no denying that FTAs and disruptions are more powerful than elaborate, well-crafted speeches.

My point here is that the use of interruptions, disruptions is a way to circumvent rationality thereby getting closer to the people, connecting with them in an more intimate way.

The closer to the people, the better the accountability.

Stéphane REVILLET 6 Insulting as (un)parliamentary practice in the British and Swedish Parliaments: a rhetorical approach, 2004. 7 Edelman Murray, The Symbolic Uses of Politics, 1972, p.31

"a piece of theatre", in Punch and Judy Politics, Ayesha Hazarika and Tom Hamilton. "A major performance at theatre", in Prime Minister's Questions as Political Ritual at Westminster, J. Lovenduski.