Quasi-Bayes properties of a recursive procedure for mixtures
Sandra Fortini Sonia Petrone



[
    [
    December 30, 2023
=====================




Bayesian methods are attractive and often optimal, yet nowadays pressure for fast computations, especially with streaming data and online learning, brings  renewed interest in faster, although possibly sub-optimal, solutions. To what extent these algorithms may approximate a Bayesian solution is a problem of interest, not always solved. On this background, in this paper we revisit a sequential procedure proposed by <cit.> for unsupervised learning and classification in finite mixtures, and developed by  M. Newton et al. (<cit.>), for nonparametric mixtures. Newton's algorithm is simple and fast, and theoretically intriguing. Although originally proposed as an approximation of the Bayesian solution, its quasi-Bayes properties remain unclear.  We  propose a novel methodological approach. We regard the algorithm as a probabilistic learning rule, that implicitly defines an underlying probabilistic model; and we find this model. We can then prove that it is, asymptotically, a Bayesian, exchangeable mixture model. Moreover, while the algorithm only offers a point estimate, our approach allows us to obtain an asymptotic posterior distribution and asymptotic credible intervals for the mixing distribution. Our results also provide practical hints for tuning the algorithm and obtaining desirable properties, as we illustrate in a simulation study. 
Beyond mixture models, our study suggests a theoretical framework that may be of interest for recursive quasi-Bayes methods in other settings.




Keywords. Asymptotic exchangeability. Bayesian nonparametrics. Conditionally identically distributed sequences. Dirichlet process. Predictive distributions. Recursive learning. 




§ INTRODUCTION

Bayesian methods have always been attractive, for their internal coherence, their rigorous way of quantifying uncertainty through probability, their optimal properties in many problems. Analytic difficulties have been overcome by  efficient computational methods and Bayesian procedures are nowadays widely and successfully used in many fields. However, fast computations remain a challenge, that hampers an even wider application of Bayesian methods among practitioners; the more so with streaming data and online learning, where inference and prediction have to be continuously updated as new data become available. 
In the modern trade off between information and computational efficiency, slightly misspecified but computationally more tractable methods receive renewed interest, as a reasonable compromise. Popular algorithms, such as Approximate Bayesian Computation (ABC) and variational Bayes, arise as approximations of an optimal Bayesian solution. Indeed, one could expect 
 that a method which performs well is at least approximately Bayes. For a Bayesian statistician, the capacity of a learning scheme to be, at least approximately, a Bayesian learning scheme should be a minimal requirement for its validation.

This is the underlying theme of this paper, that we study and illustrate in a specific, important, case, namely sequential learning in mixture models. 
We revisit a recursive algorithm initially proposed by <cit.> for unsupervised sequential learning and classification in finite mixtures, and extended by M. Newton and collaborators  
(<cit.>, <cit.>,  <cit.>, <cit.>), 
to provide a recursive fast approximation of the computationally expensive Bayesian estimate of the mixing density in nonparametric mixture models.
Recent interesting developments are in <cit.>. Convergence results have validated these algorithms as  consistent frequentist estimators (<cit.>, <cit.>, <cit.>, <cit.>). However, to what extend they  provide an approximation of a Bayesian procedure, as for the original motivation, is not fully understood. In this paper we address this question, shedding light on their  quasi-Bayes properties.  Beyond mixture models, we believe that our results may offer a methodological viewpoint of interest for quasi-Bayes recursive computations in other settings. 

Let us start by reminding the motivating problem in <cit.>, sequential unsupervised learning and classification by mixtures. 
One wants to recursively classify observations x_1, x_2, … in one of k populations (e.g., pattern types, or  signal sources, etc.), 

with no feedback about correctness of previous classifications. A finite mixture model for this task assumes 

    X_n |∑_j=1^k π_j   f_j(x).

Here the mixture components f_j(·) are known (extensive studies may be available on the specific components), but the mixing proportions =(π_1, …, π_k) are unknown. The classical Bayesian solution assigns a Dirichlet prior distribution on the unknown proportions, ∼(α_1, …, α_k), and proceeds by Bayes rule. 
Learning is solved through the posterior distribution p(| x_1, …, x_n) and classification through the predictive probabilities that X_n+1∼ f_j, j=1, …, K,  given (x_1, …, x_n). 
Unfortunately, computations, especially sequential computations, are involved.

The above model is a special case, with a discrete mixing distribution  having  atoms 1,…, k and unknown masses π_1, …, π_k, of a general mixture model 

    X_i | f_(x)=∫ f(x |θ) d (θ).

Here, a problem of interest is recursive estimation of the mixing distribution .  
The mixture model (<ref>) can be equivalently expressed in terms of a latent exchangeable sequence (θ_i), such that the X_i are conditionally independent, given (θ_i), with X_i |θ_i ∼ f(x |θ_i), and the θ_i are a random sample from , that is, θ_i |. Then, the Bayesian estimate of  with respect to quadratic loss, E((·) | x_1, …, x_n), coincides with the predictive distribution P(θ_n+1∈·| x_1, …, x_n) of θ_n+1 given the data. Thus, the problem can be rephrased as recursive prediction. 

In a Bayesian nonparametric approach, a prior with large support is assigned on the random mixing distribution ; typically, a Dirichlet process (DP), with parameters α and G_0, ∼(α, G_0). mixtures have expanded in  a floury of applications in many areas, and although 
many extensions have been developed, they remain a basic reference in Bayesian nonparametrics.  
Computational difficulties can be addressed by MCMC methods, or by   variational Bayes or ABC approximations. 
If the observations x_i arrive sequentially, one may resort to sequential Monte Carlo methods, sequential importance sampling (<cit.>), or more recent sequential variational Bayes methods (<cit.>, <cit.>), or combination of them (<cit.>). 
Still, these methods have a computational cost (for example, in the optimization steps), or only have a heuristic derivation;  simple and fast recursive algorithms remain attractive. 

Motivated by computations in mixture models, M. Newton and collaborators 
(<cit.>; <cit.>, <cit.>, <cit.>) proposed a simple recursive  estimate for , that starts at an initial guess G_0 and for n ≥ 1 recursively computes
the estimate 

    G_n(A) = (1-α_n) G_n-1(A) + α_n 
    ∫_A f(x_n |θ_n) d G_n-1(θ_n)/∫_Θ f(x_n |θ_n) dG_n-1(θ_n),

where (α_n) is a sequence of real numbers in (0,1) and it is usually  assumed that α_n→ 0 as n→∞, with ∑_nα_n=∞ and ∑_nα_n^2<∞. A standard choice, in analogy with DP mixtures, is α_n=1/(α+n) with α>0.  
For finite mixtures, as in (<ref>),  the rule (<ref>) corresponds to the sequential procedure of Smith and Makov (1978).  

<cit.> propose the recursive rule (<ref>) in applications to interval censored data and mixtures of Markov chains, further developed and extended, together with theoretical properties, by <cit.>. Theoretical properties have been studied from a frequentist viewpoint: that is, regarding G_n as an estimator of the mixing distribution, under the assumption that the data are  independent and identically distributed (i.i.d.) according to a true (identifiable) mixture model f_G_true.  <cit.> prove frequentist consistency of their recursive 
estimator, for finite mixtures, using stochastic approximation techniques. 
<cit.> shed light on the connection with stochastic approximation, thus relating frequentist consistency of the algorithm to the convergence properties of stochastic approximation sequences. 
They prove consistency of G_n for a discrete G_true with known atoms, and extend to the case of mixture kernels with an unknown common parameter. <cit.> and <cit.> prove frequentist weak consistency of the estimator (<ref>), under conditions on the mixture kernels, and give results on convergence in probability for a permutation-invariant version of it. 

Again, these results regard Newton's algorithm (<ref>) as a frequentist estimator. Its theoretical properties as an approximation of the computationally expensive Bayesian solution remain quite unsolved.  
One could argue that, when consistent, Newtons' estimator will asymptotically agree, almost surely with respect to the law of i.i.d. observations from a true F_G_true  (F_G_true^∞-a.s.), with a consistent Bayesian estimator for . 
But, also, with any other consistent estimator for ; which flattens their different nature. 
These results are important, and Newton's recursive estimator has the advantage of being  computationally faster than other consistent estimators for the mixing distribution; but its Bayesian motivation is lost. 
 
We take another, we believe more properly Bayesian, approach. Our starting point is that, when using (<ref>), a researcher is changing the learning rule on θ_n+1;  therefore, implicitly, using a  probabilistic model that is different from the Bayesian exchangeable model (<ref>). What is this model? Is it quasi-Bayes? This remark potentially applies to other approximation algorithms  
which, as well, more or less implicitly, use a probabilistic model different from the stated Bayesian one.

To address these questions, we need to formalize a notion of a quasi-Bayes procedure.  
The term quasi-Bayes is used under many meanings in the literature (see e.g. <cit.>). We mutuate this term from <cit.>, and formalize its meaning as follows. 
Details are given in Section <ref>. Mixture models specify an exchangeable probability law for the infinite sequence (X_n). 
A quasi-Bayes learning rule for mixtures should preserve this invariance property, at least asymptotically. This holds more generally, beyond mixture models. A first requirement we ask for a probabilistic learning rule to be a quasi-Bayes approximation of a Bayesian exchangeable learning rule is asymptotic exchangeability. Informally, the observations' labels do not matter if we look at far enough pieces of the sequence (X_n). 
Yet, asymptotic exchangeability is only a minimal property. A refinement refers to a specific exchangeable law. 
Let   be an exchangeable probability law for (X_n). We say that P is a quasi-Bayes approximation of    if it is asymptotically exchangeable,  and the exchangeable limit sequence has probability law .  

On this basis,  we address the following questions: 

1. If  a researcher uses the recursive procedure (<ref>) as a probabilistic learning rule, that is, she uses (<ref>) as the predictive distribution for the latent parameter θ_n+1, given x_1:n, what probabilistic model is she implicitly assuming for the observable (X_n)?  Is it an approximation, at least asymptotically, of a Bayesian, exchangeable, mixture model? 

We prove that the probabilistic model underlying the recursive 
rule (<ref>) is indeed asymptotically exchangeable. More precisely, 
it implies that the sequence (X_n) is conditionally identically distributed (c.i.d.) (<cit.>, <cit.>); roughly speaking, for any n, future observations X_n+k, k ≥ 1, are identically distributed, given X_1:n. 
 For stationary sequences, the c.i.d. property is equivalent to exchangeability; in general, a c.i.d. sequence is not exchangeable, but it is such asymptotically. This first result says that a researcher using (<ref>) as the predictive rule is implicitly assuming some form of non stationarity in the data, that tends to vanish  in the long run. A c.i.d. model could actually be the appropriate model in situations where exchangeability is broken by competition, selection or other forms of non stationarity, but the system converges to a stationary, exchangeable steady state. If, instead, it is used as an approximation of a honest exchangeable model, it guarantees the minimal property of being asymptotically exchangeable: (X_n+1, X_n+2, …) converge in distribution to an exchangeable sequence, say (Z_n). See Section <ref>.

Asymptotic exchangeability alone does not explicitly provide the statistical model of the limit sequence (Z_n); which is also, by the above result, the asymptotic statistical model for the c.i.d. sequence (X_n). We thus refine the  result by finding such model. Namely, we prove that there exist a random distribution G, which is the almost sure weak limit of the sequence G_n, such that the Z_n are conditionally i.i.d. given G, according to a distribution F_G with density f_G of the form (<ref>). Therefore, so are the X_n, asymptotically; roughly speaking, X_i | G  F_G for n large, where  means approximately i.i.d.
In  this sense, Newton's recursive learning rule arises from a quasi-Bayes mixture model.

These results shed light on an open question by <cit.>. Although Newton's recursive rule is motivated by computations in  mixture models, they show two examples where the Bayesian estimate of G with a prior and the recursive estimate G_n have different performance. 
Thus, they pose an open question: If Newton's recursive algorithm is not an approximation of the prior Bayes estimate, for what prior does the recursive estimate approximate the corresponding Bayes estimate? Our results explain the reason why the random distribution G may have a probability law  far apart from a DP. The recursive rule (<ref>) underlies a sequence (θ_n) which is not exchangeable; thus, there is no random mixing distribution G such that θ_i | G  G (one might rather think, in a state-space fashion, of a sequence of random distributions _n, such that  θ_n |_n _n; see Section <ref>). However, we show that such G exists asymptotically, 
and θ_n | G  G for n large. The probability law of G (the prior”) is implicitly determined by the c.i.d. sequence (X_n), through its so-called directing random measure F_G; and, if the mixture F_G is identifiable, it is unique. (Notice that we are denoting with the same symbol a distribution function (d.f.)  and the corresponding measure). 
Results on the explicit law of the directing random  measure of c.i.d. sequences are known only in limited, simple cases. We do obtain 
that the directing random  measure for (X_n) is a mixture of the form F_G, but we do not have explicit results on the distribution of G. However, we can prove that, under fairly mild conditions, the random distribution G is absolutely continuous, a.s.; thus, its law is rather a smoothing of a Dirichlet process, the latter a.s. selecting discrete distributions. In fact, Newton's algorithm was originally given in terms of densities, assuming that the unknown mixing distribution is absolutely continuous, with density g. 
Notice that, if G_0 is absolutely continuous with respect to a sigma-finite measure λ on Θ  (denoted G_0≪λ), 
then also G_n≪λ, a.s., for every n ≥ 1, and its density g_n 
satisfies the recursive equations

    g_n(θ)   =(1-α_n)g_n-1(θ)+α_ng_n-1(θ)f(x_n|θ)/∫ f(x_n|θ')g_n-1(θ')dλ(θ').

We have given Newton's rule (<ref>) for the distributions, as they are better defined objects when we read the algorithm as a probabilistic learning rule.   

2. The second main question we address is as follows. As an algorithm, Newton's recursive rule (<ref>) only gives a point estimate of the mixing distribution G. Can one provide a more complete description of the uncertainty? 
A Bayesian approach would fully describe the uncertainty through the posterior distribution. Can one enrich Newton's algorithm, by providing a posterior distribution on G?

Our key to address this question is, again, to regard Newton's rule (<ref>) as a probabilistic learning rule. We can then formally prove that G_n(·)=E(G(·) | x_1, …, x_n), 

thus, it is, indeed, the point estimate of G under quadratic loss; moreover, although the prior distribution of G is only implicitly defined, we can approximate the posterior distribution of G, that results from such implicit prior. More specifically, we provide convergence rates and an asymptotic Gaussian approximation of the finite dimensional conditional distribution P([G(A_1), …, G(A_k)] ∈·| x_1, …, x_n), for any fixed measurable sets A_1, …, A_k, k ≥ 1. 
Thus, not only one has a quasi-Bayes point estimate, but may provide asymptotic credible regions. 



These results illuminate on the probabilistic model implied by Newton et al.'s recursive rule. Interestingly, this understanding also gives new insights on the role of the weights α_n and of other main parameters of the model. As we illustrate in Section 5, it provides practical advise about tuning the algorithm for obtaining desirable properties; notably, for attenuating the effect of the lack of exchangeability on the estimates, a problem usually addressed, with higher computational cost, by taking averages over permutations of the sample. 



Section <ref> recalls the basic structure of DP mixture models, and introduces our viewpoint on Newton et al.'s recursive algorithm as a probabilistic learning rule. We prove its quasi-Bayes properties in Section <ref>. Section <ref> provides rates and asymptotic results, together with asymptotic credible intervals. A simulation study for location mixtures of Gaussian distributions, in Section <ref>, illustrates how the undesirable sensitivity of the estimates to the ordering of the observation can be attenuated, by tuning the weights and the main parameters of the model.  All the proofs are collected in Section <ref>. Directions for developments are finally discussed.




§ DIRICHLET PROCESS MIXTURES AND A NEW LOOK AT NEWTON'S SEQUENTIAL PROCEDURE



We first recall the basic structure of Bayesian inference for DP mixture models, in order to motivate in more details the recursive rule (<ref>), and to introduce some further notation. 

Again, the DP mixture model has a 
hierarchical formulation in terms of a latent exchangeable sequence (θ_i)

    X_i |θ_i  f(x |θ_i)  
          θ_i |, 
          ∼(α G_0) ,
   
where (<ref>) is a short notation for 
X_n | X_1,…,X_n-1, (θ_n) ∼ f(x|θ_n), for every n ≥ 1, and f(·|θ) is a density with respect to a sigma-finite measure μ on the sample space 𝕏. Integrating the θ_i out, one has the mixture model (<ref>), with a  prior on .  We denote by  the probability law on the process ( (X_n, θ_n)) so defined.  
We assume that the X_i and the θ_i are, respectively, random variables with values in 𝕏⊆ℝ^d and Θ⊆ℝ^p, equipped by their Borel sigma-fields ℬ(𝕏) and ℬ(Θ);  but the results hold more generally, for Polish spaces. 
Throughout the paper, we refer to conditional distributions as regular versions. We use the short notation X_1:n=(X_1, …, X_n), and P(A | x_1:n) for 
P(A | X_1=x_1, …, X_n=x_n). 
Unless explicitly stated, convergence of distributions is in the topology of weak convergence. Weak convergence of F_n to F  is denoted by F_n⇒ F.

Inference on  in a mixture model moves from the conditional distribution  

    |θ_1:n, x_1:n∼(α G_0 + ∑_i=1^n δ_θ_i),
 
to get the posterior distribution as a mixture of DPs

    | x_1:n∼∫(α G_0 + ∑_i=1^n δ_θ_i)    d (θ_1:n| x_1:n).

The Bayesian point estimate _n^(Bayes) of , with respect to quadratic loss,  is the conditional expectation of , and coincides with the predictive distribution of θ_n+1, given X_1:n. 
By the urn structure characterizing the Dirichlet process

    (θ_n+1∈·|θ_1:n, x_1:n) = α G_0(·) + ∑_i=1^n δ_θ_i(·)/α+n,

therefore

    _n^(Bayes)(·)   =   E((·) | x_1:n)= (θ_n+1∈·| x_1:n)     
       =   α G_0(·) + ∑_i=1^n (θ_i ∈·| x_1:n)/α+n
       =   α+n-1/α+n α G_0(·) + ∑_i=1^n-1(θ_i ∈·| x_1:n)/α+n-1 + 1/α+n__n-1^Bayes(θ_n ∈·| x_n),

where  

    __n-1^(Bayes)(θ_n ∈ A | x_n) = 
    ∫_A f(x_n |θ) d _n-1^(Bayes)(θ)/∫_Θ f(x_n |θ) d _n-1^(Bayes)(θ).

In the Bayesian estimate, as a new observation x_n becomes available, the information on all the past θ_i, i=1,…, n-1, is updated. This efficiently exploits the sample information, but is computationally expensive. Instead, Newton's algorithm (<ref>) does not update the estimate G_n-1, and x_n only enters in inference on θ_n, in an empirical Bayes flavor. The two estimates coincide only for n=1. However, even for n=1, Newton's rule makes a simplification of the posterior distribution of G, replacing the mixture of Dirichlet processes  ∫(α G_0+δ_θ_1) d(θ_1 | x_1), as from  (<ref>),  with a (α G_0 + (θ_1 | x_1)). For n ≥ 1, Newton's estimate loses in efficiency, not fully exploiting the sample information, but, on the other hand, is very fast; if one evaluates (<ref>) on a grid of m points and calculates the integral that appears at the denominator using, say, a trapezoid rule, then the computational complexity is m n. 
 


 §.§ Newton's algorithm as a probabilistic learning rule


As underlined, in DP mixture models, the point estimate E(| x_1:n) corresponds to the predictive distribution of θ_n+1| x_1:n. Our point is that, similarly, Newton's recursive rule should  be regarded as a probabilistic learning rule for model (<ref>), that expresses a different, computationally simpler, predictive distribution for the θ_n: 

    [             X_n |θ_n  f(x |θ_n); θ_n+1| x_1:n∼ G_n(·),    n ≥ 1, ]
   
with θ_1 ∼ G_0 and G_n given by (<ref>). 

A known result in probability theory is that the predictive rule for a sequence of random variables characterizes the probability law of the sequence (Ionescu-Tulcea Theorem). On this basis, a researcher using the predictive rule (<ref>) is implicitly using a different probabilistic model for the sequence ((X_n, θ_n)), in place of the exchangeable mixture model (<ref>). Our aim is to make this model explicit. Such model may be of autonomous interest in some experimental circumstances.
When we regard Newton's sequential procedure as a probabilistic learning rule, it provides the formal framework for understanding its quasi-Bayes properties. 

Let us denote by P a probability law on ((X_n, θ_n)) that is consistent with assumptions (<ref>). Newton's recursive formulae can now be read as probabilistic implications, under the law P. The estimate G_n can be written in a prediction-error correction form 

    G_n(·)=G_n-1(·)+α_n [P(θ_n∈·| x_1:n)-P(θ_n∈·| x_1:n-1)],

where the correction term can now be interpreted as a difference between predictive distributions, computed according to P. Moreover, we can appreciate the different information conveyed in the predictive rule, with respect to mixtures. Simple computations show that one can write G_n as  

    G_n(·)=α G_0(·)+∑_k=1^nγ_k P(θ_k∈·| x_1:k)/α+∑_k=1^nγ_k,      n≥ 1,

where P(θ_k∈·| x_1:k)= P_G_k-1(·| x_k);  α>0, γ_1=α_1α/(1-α_1) and γ_n=α_n(α+∑_k=1^n-1γ_k)/(1-α_n) for n≥ 2. For α_n = 1/(α+n), one has γ_n=1 for all n ≥ 1.   In this case, one has a direct comparison with the corresponding formula (<ref>) for DP mixtures.  The latter originates from the urn structure characterizing the Dirichlet process. This suggests that Newton's recursions are based on a different urn scheme, possibly urns of distributions; see <cit.>. 

 
A further, immediate implication of the assumptions (<ref>) is that, under the law P, the θ_i, and consequently the X_i, are no longer exchangeable. In fact, we show that Newton's model (<ref>) replaces exchangeability with  a weaker form of dependence; namely, with (X_n) 
being a conditionally identically distributed sequence. 
This is an interesting stochastic dependence structure, noticed by <cit.>, and developed by <cit.>. 
It implies that the  X_n are  asymptotically exchangeable. 
Before proceeding, let us remind some basic definitions and properties of c.i.d. sequences. 



 §.§ Conditionally identically distributed sequences

 
<cit.> (Proposition 2.1) 
proves that a stationary sequences that satisfies 

    (X_1, …, X_n, X_n+2) d=(X_1,…,X_n,X_n+1),    n ≥ 1

is exchangeable. Clearly, the converse is true, thus condition (<ref>) is equivalent to exchangeability for stationary sequences. 
Notice that (<ref>) implies that (X_1, …, X_n, X_n+k) d=(X_1,…,X_n,X_n+1), for any n ≥ 1 and k≥ 1; informally, for any n ≥ 1, 

    X_n+k| X_1:n X_n+1| X_1:n,   ,
 
where  means equal in distribution.  <cit.> extend this notion  and introduce the term conditionally identically distributed sequences with respect to a filtration.
 
Let = (_n)_n ≥ 1 be a filtration. A sequence of random variables (X_n) is conditionally identically distributed with respect to the filtration  (-c.i.d.) if it is adapted to  and, a.s., 

    E[ h(X_n+k) |ℱ_n] = E[ h(X_n+1) |ℱ_n],

for all k ≥ 1, n ≥ 0 and all bounded measurable functions h: 𝕏→ℝ 
 
Less formally, a sequence (X_n) adapted to  is -c.i.d. if the X_i are identically distributed and X_n+k|_n  X_n+1|_n,    k ≥ 1, n ≥ 1.  
When  is the natural filtration of (X_n), the sequence is said to be c.i.d. An -c.i.d. sequence is also c.i.d.  
C.i.d. sequences preserve main properties of exchangeable sequences. In particular, as for exchangeable sequences, the sequence of the empirical distributions and the sequence of the predictive distributions  converge a.s., to the same random distribution, F say. If (X_n) is c.i.d. with probability law P, then 

    F̂_n ≡∑_i=1^n δ_X_i/n⇒ F     
    P_n ≡ P(X_n+1∈·| X_1:n) ⇒ F ,   .
 
For exchangeable sequences, the limit F is called the directing random measure  (the statistical model, in Bayesian inference) and the probability law of F is the de Finetti measure (the prior distribution). The limit F in (<ref>) is referred as the directing random measure for c.i.d. sequences, too.

An exchangeable sequence is clearly c.i.d., but the reverse is not generally true. However, c.i.d. sequences are asymptotically exchangeable. 

  
A sequence of random variables (X_n) is asymptotically exchangeable, with directing measure F, if  

    (X_n+1, X_n+2, …) d→ (Z_1, Z_2, …)
  
for an exchangeable sequence (Z_n), with directing random  measure F.  


For a sequence (X_n), convergence of the predictive distributions to a random probability measure, μ say,  implies that the sequence is asymtotically exchangeable, with directing random measure μ
(<cit.> Lemma 8.2).  
Thus, by (<ref>), a c.i.d.  sequence (X_n) is asymptotically exchangeable, with directing random measure F. Informally, X_n | F ind≈ F, for large n. 

Fundamental results and limit theorems for -c.i.d. sequences are given by <cit.>. Applications in Bayesian nonparametric inference include <cit.> and the c.i.d. hierarchical model proposed by <cit.>. 



§ QUASI-BAYES PROPERTIES OF NEWTON'S MODEL
 

Newton's model (<ref>) does not fully specify the probability law of the process ( (X_n, θ_n) ), as it only assigns the predictive distribution of θ_n+1 conditionally on the observable X_1:n, while not enough restrictions are made on the conditional distribution of θ_n+1 given X_1:n and θ_1:n. Still, model (<ref>) has interesting implications, that we study in this section. Clearly,  an obvious way to obtain a full specification is to additionally assume that θ_n+1 is conditionally independent on  θ_1:n, given X_1:n. This stronger assumption  might be motivated by  the non-stationary nature of the sequence (θ_n) in (<ref>), and would considerably simplify the analysis; yet, it is not necessary, and our results are developed under the only assumptions  (<ref>). 



The first remark is that Newton's model (<ref>) implies that the (X_n) are no longer exchangeable. We will show that (X_n) is c.i.d. In some applications, exchangeability may actually be broken by forms of non-stationarity, and a c.i.d. model may offer a sensible description of the phenomena.  

Example. As a simple example, suppose that data are observed over time, after an intervention that affects the population under investigation, but whose effect is unpredictable' and tends to vanish. One still assumes that X_n |θ_n  f(x_n |θ_n), but a disequilibrium is introduced, so that the sequence (θ_n) is no longer exchangeable. 
One may rather envisage, in a state-space fashion, a sequence of random distributions, (_n) say, such that θ_n |_n _n. Unpredictability of the dynamics may be expressed by further assuming that, for any n ≥ 1,  the conditional law of the new random distribution _n, given θ_1:n-1, x_1:n-1, only depends on x_1:n-1 and is Dirichlet process, centered on the current estimate G_n-1  

    _n  | X_1:n-1, θ_1:n-1∼( (α+n-1) G_n-1).

Then, as x_n becomes available, Newton's one-step-ahead update (<ref>) is exact, that is, it is the Bayesian estimate of _n from the prior (<ref>) 

    P(θ_n+1∈·|θ_1:n, x_1:n) =
    E (P(θ_n+1∈·|_n) |θ_1:n, x_1:n) 
    = E( _n |  x_1:n) = G_n .
 
As we will show, the sequence G_n converges a.s. to a random distribution G. By properties of the DP, this implies that the conditional law (<ref>) of _n converges a.s. to a measure degenerate on G. This fact expresses the idea that the disequilibrium tends to vanish. Intuitively,  θ_n | G  G for n large. We will indeed prove that the θ_n are asymptotically  exchangeable. 



In more standard setting, stationarity holds; the  researcher would regard the X_i as exchangeable, judging that the order of the observations does not matter. The lack of exchangeability implied by assumptions (<ref>) is thus a misspecification, only motivated by the need of fast computations. Then a minimal requirement is that (X_n) is at least asymptotically exchangeable; informally, the order does not matter if we look at X_n, X_n+1, …, for n large enough. Moreover, one would want a mixture model of the form X_n | G  f_G, at least  asymptotically. To prove the latter property, we start by showing that such G exists, and is the a.s. limit of the sequence of the predictive distributions G_n. Furthermore, as for exchangeable sequences, G_n is the expectation, given X_1:n, of such limit distribution G.
The proof of the theorem below, as well as the proofs of all the following results, are collected in Section <ref>.  

 
Let the process ((X_n, θ_n)) have probability law P that satisfies assumptions (<ref>). Then 


  (i)	the sequence (G_n) converges to a random probability measure G, a.s.; 

  (ii) 
for every n ≥ 1 and measurable set A, 
P(θ_n+k∈ A | X_1:n)= E( G(A) | X_1:n), for all k ≥ 1. 


 
An immediate consequence of the weak convergence of G_n to G is that ∫ h(θ)dG_n(θ) →∫ h(θ) dG(θ) a.s. for any continuous and bounded function h on Θ. We prove  that the convergence can be extended to functions h that are  integrable with respect to G. 

	
Let ((X_n, θ_n)) ∼ P satisfy the  assumptions (<ref>), and let h(·) be a measurable function on Θ, such that 
∫ |h(θ)| dG(θ) < ∞ a.s. 
Then, for n →∞, 

    ∫ h(θ) dG_n(θ) 
    →∫ h(θ) dG(θ),

The condition ∫ |h(θ)| dG(θ) < ∞ a.s. holds, in particular, if h is measurable and ∫ |h(θ)| dG_0(θ) < ∞.

 
The following theorem proves that Newton's learning rule (<ref>) implies that the sequence (X_n) is c.i.d., thus asymptotically exchangeable, and that its directing random measure has a mixture density of the form f_G. 
In this sense, Newton's model is a quasi-Bayes mixture model. 

  	
Let ((X_n, θ_n)) ∼ P satisfy  the assumptions (<ref>). 
	Then 


  (i) The sequence (X_n) is c.i.d.; 

  (ii)  The sequence of predictive densities f_G_n converges in L_1 to 
f_G ≡∫ f(x |θ) dG(θ), P-a.s.,  where G is the a.s. weak limit of (G_n); 
 
  (iii) 
(X_n) is asymptotically exchangeable, and its directing random measure has density f_G with respect to μ.


 
Informally, the above results 
say that X_n | G  f_G, for n large. Notice that G plays the role of the (infinite-dimensional) parameter of the directing random measure F_G of (X_n), and as such, it is a function of (X_n). If the mixture is identifiable, F_G uniquely determines G. 
Moreover, by properties of c.i.d. sequences, F_G is also the a.s. weak limit of the sequence of empirical distributions ∑_i=1^n δ_X_i/n. Therefore, it is measurable with respect to the tail sigma-field of (X_n). 

Intuitively, asymptotic exchangeability of the sequence (X_n) implies that  (θ_n) is also asymptotically exchangeable. The following theorem provides a formal proof. In fact, if we assume the additional condition that θ_n+1 is independent on θ_1:n, given  X_1:n, then (θ_n) is easily proved to be c.i.d., thus asymptotically exchangeable. Again, the theorem below does not use this assumption. 
 

	
If the mixture F_G=∫ f(y |θ) dG(θ) is identifiable, then Newton's learning scheme (<ref>) implies that the sequence (θ_n) is asymptotically exchangeable, with directing random measure G corresponding to the a.s. limit of the sequence G_n.

 



 §.§ On the prior distribution of G


In its original derivation, Newton's algorithm, although moving from mixtures, is proposed as a recursive estimate of the mixing density, implicitly assuming that the unknown  mixing distribution G is absolutely continuous, with density g. Regarding Newton's algorithm as a probabilistic learning rule, we have shown that it implies a quasi-Bayes mixture model, where, asymptotically, the X_i are a random sample from a mixture f_G. Yet, compared with mixtures, Newton's model changes the prior on G, which is no longer (in general) a DP. 
Giving explicit results on the law of the limit G is difficult. In the literature on c.i.d. processes, there are very few results, for very simple cases. Even for exchangeable Bayesian nonparametric methods, giving explicit results on the prior implied by predictive constructions often is a non trivial problem.  
Although we cannot explicit the prior distribution on G implied by Newton's model (<ref>), we can prove that, under fairly mild sufficient conditions, G is absolutely continuous, a.s. Moreover, although the prior is only implicitly defined, in the next section we give an asymptotic Gaussian approximation of the (finite-dimensional) posterior distribution that results from such prior.   

As noticed, if G_0 is absolutely continuous with respect to a sigma-finite measure λ on Θ,  G_0 ≪λ, then G_n ≪λ, and its density g_n satisfies Newton's recursive rule (<ref>). It is immediate to verify that, for any fixed θ,  the sequence (g_n(θ)) is a martingale.  Since g_n(θ) is non-negative, there exists a function g^*(θ) such that, for every θ,  g_n(θ) converges to g^*(θ) a.s. However, this fact is not sufficient to conclude  that G ≪λ. The latter property requires that g_n converges in L^1, or, equivalently, G_n converges to G in total variation. This fact is made precise in the following Lemma <ref>. 
Then, Lemma <ref> gives sufficient conditions for G_n → G in total variation, and G≪λ, a.s.
The two lemmas are extensions of Theorems 1 and 4 in  <cit.>, the main difference being that we do not assume that the involved sequence of random variables is c.i.d. To be -c.i.d.,  a sequence has to be adapted to the filtration , and this condition is not satisfied, in general, by the sequence (θ_n) in Newton's model.  
Still, one can fairly easily extend the proofs of Theorems 1 and 4 in <cit.>, by directly requiring the martingale property of the sequence of random measures Q_n, which is otherwise implied by the -c.i.d. property, to obtain the two lemmas below, whose proofs are therefore omitted. The proofs of the following Theorem <ref> is instead provided in Section <ref>.


	
	Let λ be a sigma-finite measure on a Polish space S. For any n, let Q_n be a random measure on S such that the sequence (Q_n) is a measure-valued martingale with respect to a filtration (ℱ_n), and let Q be its limit. 
	Then Q ≪λ a.s.  if and only if, a.s.,   
Q_n ≪λ for every n and Q_n converges to Q in total variation.
	


	
	Let λ, Q_n and Q be as in Lemma <ref>.
Assume that Q_n ≪λ a.s. for every n, with density  q_n. 
Then Q≪λ a.s.  if and only if for every compact K such that λ(K)<∞, 

q_n is, a.s., a  function on S uniformly integrable with respect to λ_K, 

where λ_K(·)=λ(·∩ K) is the restriction of λ on K. 

In particular, Q≪λ a.s. if, for every K compact, there exists p>1 such that, a.s., 

    sup_n∫_K q_n(x)^pdλ(x)<∞.

A sufficient condition for (<ref>) is 

    sup_nE(∫_K q_n(x)^pdλ(x))<∞.


  
  
  
We can now provide sufficient conditions for the a.s. absolutely continuity of the limit mixing distribution G in Newton's model.


	
Let G be the a.s. limit of the sequence of predictive rules G_n defined by (<ref>). If the following conditions hold

    ∑_nα_n^2<∞  ;      ∫_Kg_0(θ)^2dλ(θ)<∞,  ;

and 

    sup_θ_1,θ_2∈ K∫f(x|θ_1)^2/f(x|θ_2)dμ(x)<∞ ,     ,

 then G≪λ, a.s.
Moreover, a.s, g_n converges in L^1 to g ≡ dG/dλ.
	
	
Assumptions (<ref>)  are quite natural. They hold, for example, if α_n=1/(α+n) and g_0 is continuous or bounded. Assumption (<ref>) is more delicate. It holds, for example, if f is a Poisson density or a Gaussian density with fixed variance or a Gamma density with fixed shape parameter. 
 A similar assumption is considered in <cit.>. 

In the next subsection, we give a simple comparison between the and the law of G arising from Newton's recursions, through a simulation study. 



 §.§ Empirical study
 
mixtures are widely used for clustering, based on the random partition of the θ_i induced by the DP. For model (<ref>), the asymptotic mixing distribution G is absolutely continuous under fairly mild conditions; thus, one rather has multiple shrinkage effects in the estimation of θ_1:n, governed by the modes of the mixing density g. 
Nevertheless, one may expect that, if x_i provides strong information on θ_i, then Newton's predictive rule G_n implies a negligible loss of information with respect to the predictive rule G_n^(Bayes) corresponding to a mixture (compare (<ref>) and (<ref>), with γ_k ≡ 1). Consequently,  in this case, the law of G should be close to the DP. 



This behavior is illustrated in Figure <ref>. We consider a  location mixture of Gaussian distributions (θ_i, σ^2), with σ^2 known.  In this case, the assumptions of Theorem <ref> are satisfied, and  G is absolutely continuous, a.s.
The prior guess G_0 is (1,3). The weights in Newton's recursions are α_n=1/(α+n), with α=5. 
Were G ∼(α, G_0), then, for any set A, the law of G(A) would be a Beta distribution, with  parameters  (α G_0(A), α (1-G_0(A)). 
The panels in the first column of Figure <ref> compare such Beta density (gray dotted line), for A=(-∞, 0], with a Monte Carlo approximation of the law of G(0) arising from the c.i.d. model. The Monte Carlo sample is obtained by generating samples x_1:N^(m), m=1, …, 200 from the c.i.d. model (<ref>), with N=1000. For each sample, we compute G_N(0), that we take as a fairly reasonable proxy of the realization of the random G(0).  The first raw in Figure <ref> corresponds to σ^2=0.01. In the second raw, σ^2=1. Clearly, the results do not contradict the conjecture that, if σ^2 is (very) small, the law of G is close to  the DP.   
The shrinkage effect that replaces the DP random partition of θ_1:N is illustrated in the right panels of Figure <ref>, which show the density estimate g_n in the two cases. Clearly, for σ^2=0.01, g_n closely describes a random partition of the θ_n.




§ ASYMPTOTIC POSTERIOR LAWS
 

 By part (ii) of Theorem <ref>, Newton's rule G_n can be regarded as a point estimate, with respect to  quadratic loss,  of the limit mixing distribution G, in a quasi-Bayes mixture model. Yet, a Bayesian mixture model would offer more than a point estimate; it would give a description of the uncertainty through the posterior distribution of G | x_1:n. 
The probabilistic framework we have provided for Newton's learning rule allows to go beyond point estimation, studying the conditional  distribution of G, given x_1:n.  We first obtain an asymptotic Gaussian approximation of the conditional  distribution of G(A), given x_1:n, for any measurable set A. We then extend the results to the joint conditional distribution of a random vector [G(A_1), …, G(A_k)]', for any measurable A_1, … A_k. 



 §.§ Asymptotic posterior distribution and credible intervals.



Let us recall that P is a probability law for (X_n), consistent with the assumptions (<ref>). Here we  give  an asymptotic Gaussian approximation of the conditional law P(G(A) ∈·| x_1:n), for a measurable set A. 
The almost sure conditional convergence involved  is a strong form of convergence (<cit.>), that implies stable convergence  (<cit.>, <cit.>, <cit.>) and convergence in distribution of the unconditional law. 

Notice that, although having a similar flavor,  these results differ from Bernstein-von-Mises types of theorems, which are stated a.s. with respect to the law P_G_true^∞.  
The law P describes an evolutionary process and the results inform about the rate of convergence of G_n(A)(x_1:n) to the limit distribution G(A)(x_1, x_2, …), for any ω=(x_1, x_2, …) in a set of P-probability one. Berstein-von-Mises results are beyond the aims of this paper, but we will give some hints in the simulation study in Section <ref>.

Let us denote by (μ, σ^2) the Gaussian law with mean μ and variance σ^2, and by  Φ(t |μ, σ^2) its d.f. evaluated at t.  A (0,0) law is interpreted as the law degenerate at zero. The p-dimensional Gaussian law will be denoted by _p, with d.f. Φ_p.   
Without loss of generality, we can assume that f_G_0(x)≠ 0 for every x∈𝕏. This implies 

    f_G_n(x)≠ 0  .

Our first result finds a sequence  (r_n) such that the conditional distribution of √(r_n) (G(A)-G_n(A)), given X_1:n, is asymptotically a zero-mean Gaussian law, with variance

    V_A = ∫_{x:f_G(x) ≠ 0} P_G(A| x)^2 
    	dF_G(x)- G(A)^2.

We remind the notation  
P_H(A | x) =∫_A f(x|θ) dH(θ)/∫_Θ f(x|θ) dH(θ), for any d.f. H on Θ. 

Before stating the theorem, we give the following Lemma. 
Let us define, for any A ∈ℬ(Θ) and n ≥ 1, 

    V_A,n =
    ∫_𝕏 P_G_n(A| x)^2
     dF_G_n(x) - G_n(A)^2.

Notice that V_A,n can be written as V_A,n= E( (P_G_n(A | X_n+1) -  G_n(A))^2 | X_1:n), expressing the prior-to-posterior variability, given X_1:n, when G_n plays the role of the prior and P_G_n(·| x_n+1) of the posterior.  

 
For any A ∈ℬ(Θ), V_A,n converges to
V_A a.s. as n →∞.   


We can then give the main results of this section. 

 
	
	Let (α_n) satisfy  ∑_nα_n=∞ and ∑_nα_n^2<∞ and let (r_n) be a  monotone sequence of positive numbers  such that r_n∼ (∑_k > nα_k^2)^-1 as n→∞. If

    √(r_n)sup_k ≥ nα_k→ 0

	 and
	 
    ∑_k≥ 1r_k^2α_k+1^4<∞,

	then, for every A∈ℬ(Θ), 
	
    P(	√(r_n) (G(A)-G_n(A)) ≤ t | X_1:n) →Φ(t | 0,  V_A),

	with V_A as in (<ref>).
	If α_n=(α+n)^-β with 1/2 <β≤ 1 and α>0, then 
	(<ref>) holds with 
	r_n=(2β-1)n^2β-1.



Assumptions (<ref>) and (<ref>) hold for most choices of (α_n) satisfying ∑_nα_n=∞ and ∑_nα_n^2<∞. In particular, if  (α_n) is definitively decreasing, then (<ref>) is a consequence of (<ref>). A sufficient condition for (<ref>) is

    α_n=(nb_n)^-1

for a sequence b_n which is definitively non increasing. Indeed, in this case 

    lim sup_n→∞r_n-1α_n
    	=lim sup_n→∞(nb_n)^-1/∑_k≥ n(kb_k)^-2≤lim sup_n→∞(nb_n)^-1/b_n^-2∑_k≥ nk^-2≤lim sup_n→∞b_n<∞.
 
	In turn, this implies that  r_n^2 α_n+1^4<(sup_kb_k)^2α_n+1^2, for n large enough, and therefore (<ref>).
	 

 If V_A(ω)=0, then Theorem <ref> gives convergence to a degenerate distribution on zero.
From the definition of V_A, 
 it is immediate to see that V_A(ω)=0 if and only if P_G(A | X)(ω) = G(A)(ω), which happens if and only if G(A)(ω) is zero or one. 
 

In Theorem <ref>, the limit variance V_A is unknown, depending on G. By Lemma <ref>, a convergent estimator is provided by V_A,n. Replacing the random V_A with its consistent estimate, to get an asymptotic distribution that allows to compute asymptotic credible intervals for G(A) (as done through Cramér-Slutzky Theorem in standard i.i.d. settings) is here delicate, given the kind of convergence we are studying. Yet, we can prove the following 

 	
Let A ∈ℬ(Θ). Then, with (r_n) defined as in Theorem <ref> and under the same assumptions, 
for almost all ω=(x_1,x_2,…) such that V_A(ω)>0,

    P(√(r_n)G(A)-G_n(A)/√(V_A,n)≤ t |  x_1:n) →Φ(t |  0,1).

If α_n=(α+n)^-β with 1/2 <β≤ 1 and α>0, then 
(<ref>) holds with 
r_n=(2β-1)n^2β-1.




 Theorems <ref> and <ref> allow to obtain  
asymptotic credible intervals. For a fixed set A, 
 Theorem <ref> gives that, for almost all ω=(x_1, x_2, …) such that V_A(ω)>0, 

    P( G_n(A) - z_1-γ/2√(V_A,n/r_n) < G(A) < G_n(A) +  z_1-γ/2√( V_A,n/r_n)| x_1:n) ≈ 1-γ,

where z_1-γ/2 is the (1-γ/2)-quantile of the standard Gaussian distribution. If V_A(ω)=0,
then Theorem <ref> implies that the limit distribution is degenerate on zero, therefore, for any ϵ > 0

    P( G_n(A) - z_1-γ/2√(ϵ/r_n) < G(A) < G_n(A) +  z_1-γ/2√(ϵ/ r_n)| x_1:n) ≥ 1-γ,

asymptotically. It follows that, for every ϵ >0, 

    [ G_n(A) - z_1-γ/2√(max ( V_A,n, ϵ)/r_n) ; G_n(A) +  z_1-γ/2√(max ( V_A,n, ϵ)/r_n)]
 is an asymptotic credible interval for G(A), of level at least 1-γ. 




 §.§ Asymptotic joint posterior distribution and credible regions.


 
We now study the joint behavior of [G_n(A_1)-G(A_1),…,G_n(A_k)-G(A_k)]', for any fixed choice of A_1,…,A_k ∈ℬ(Θ).  
 As in the previous section, we assume that f_G_0(x)≠ 0 for every x∈𝕏, which implies (<ref>).   
For every n≥ 1, and A_i,A_i', let 

    C_A_i,A_i',n   =Cov( P(θ_n+1∈ A_i| 
    X_1:n+1),P(θ_n+1∈ A_i'|  X_1:n+1)| X_1:n)
       =∫_𝕏P_G_n(A_i | x)P_G_n(A_i'| x)dF_G_n(x)
    
    
    -G_n(A_i)G_n(A_i').
  
Furthermore, let  

    C_A_i,A_i'=∫_{x:f_G(x)≠ 0}P_G(A_i | x)P_G(A_i'| x)dF_G(x)
    -G(A_i)G(A_i').

Following the same line of reasoning as in Lemma <ref>, it can be proved that, as n→∞,
	
    C_A_i,A_i',n→C_A_i,A_i'   a.s.

Thus, the conditional covariance matrix 

    C_n(A_1,…,A_k):=[C_A_i,A_i',n]_i,i'.

converges a.s. to 

    C(A_1,…,A_k):=[C_A_i,A_i']_i,i'.

Then we can prove 

	
In model (<ref>), let (α_n) satisfy  ∑_nα_n=∞ and ∑_nα_n^2<∞.
Let (r_n) be a monotone sequence of positive numbers such that r_n∼(∑_j>nα_j^2)^-1 as n→∞.
If (<ref>) and (<ref>) hold, then,
for every k ≥ 1 and every A_1,…,A_k∈ℬ(Θ), 

    P (√(r_n)[ G(A_1)-G_n(A_1);               ⋮; G(A_k)-G_n(A_k) ]∈·| X_1:n ) →Φ_k( ·| 0,  C(A_1,…,A_k)),   ,

with C(A_1,…,A_k) as in (<ref>).
	If α_n=(α+n)^-β with 1/2 <β≤ 1 and α>0, then 
	(<ref>) holds with 
	r_n=(2β-1)n^2β-1.





The following result is the analogous of Theorem <ref> for the joint posterior distribution.


	
	Under the same assumptions as in Theorem <ref>,
for every k ≥ 1 and every A_1,…,A_k∈ℬ(Θ), 

    P(√(r_n)  C_n(A_1,…,A_k)^-1/2[ G(A_1)-G_n(A_1);               …; G(A_k)-G_n(A_k) ]∈·| x_1:n)→Φ_k(·| 0, I)   ,

for almost all ω=(x_1, x_2, …) such that (C(A_1,…,A_k)(ω)) ≠ 0. 


Based on Theorems <ref> and <ref>, we can give an asymptotic credible region for [G(A_1),…,G(A_k)]'.


 
Let  𝐆_n(𝐀)=[G_n(A_1),…,G_n(A_k)] and let χ^2_1-γ denote the (1-γ)-quantile of the chi-square distribution with k degrees of freedom. Then, for every ϵ >0, the set

    E_n^(ϵ)={𝐬∈ℝ^k:(𝐬-𝐆_n(𝐀))^T( C_n(A_1,…, A_k)+ϵ I)^-1(𝐬-𝐆_n(𝐀))
    ≤χ^2_1-γ/r_n}

satisfies

    lim inf_n P([G(A_1),…,G(A_k)]'∈ E_n^(ϵ)|
    X_1:n)≥ 1-γ   a.s.





§ SIMULATION STUDY: TUNING THE WEIGHTS Α_N
   
The lack of exchangeability of Newton et al.'s rule implies that the mixing density estimate depends  on the ordering of the observations. Indeed, the authors propose to use an average of the estimates obtained over random permutations of the original sample; see also <cit.>. Our study of the probabilistic model underlying Newton's learning rule  allows us to give new hints about how one can attenuate the effect of the ordering on the estimates. 

Figure <ref> shows the recursive estimate g_n of the mixing density, given by  (<ref>), for different permutations of the observations. 
In the first column, the data is a random sample x_1:1000, simulated from a location mixture of Gaussian distributions, with σ^2=1. 
 The true mixing distribution used in the simulation is 0.3 (-1, 2) + 0.7 (3, 1.5)
(plotted in dark gray), describing two groups in the θ_i. 
The initial distribution G_0 is (1,9), fairly diffuse (gray, dashed).  In the second and third column, the data are two different permutations of x_1:1000. 






The black lines in the plots are the estimates g_n, for n=1000. The dashed light-gray to gray  lines are the density estimates g_n for n=1; n=300 and n=500, respectively.
In the first row, the weights are α_n=1/(α+n), with α=1. The density estimates change sensitively for different  permutations of the sample. We expect that a larger value of α may attenuate the sensitivity of the estimates to permutations, for the following reason. 
We have shown that Newton's rule underlines an evolutionary, c.i.d. model, where the order (time, say) matters. Intuitively, the effect of the ordering is attenuated, if such evolution is milder. The example in Section <ref> helps the intuition. The θ_n are not exchangeable, thus, there is no random latent distribution G such that θ_n | G  G. In the example, we envisage a sequence (_n) of random distributions, that evolve over time from the initial distribution G_0. Such dynamics from G_0 is stronger, the larger is α_n; thus, when α_n=1/(α+n), the smaller is α.  This is evident by looking at (<ref>): the variance of a DP increases, the smaller is the scale parameter.  
Thus, a large value of α not only expresses that _1 ∼ DP(α G_0) is concentrated around the (absolutely continuous) initial distribution G_0 (which is reasonable, as, in fact, Newton's algorithm is thought for an absolutely continuous mixing distribution); but also, and we believe interestingly, it expresses a moderate evolution of _n; thus, a situation closer to exchangeability. The second row of Figure <ref> shows the density estimates g_n with α_n=1/(α+n) and α=100. The data are the same as in the first row. The effect of the ordering of the observations is clearly attenuated. 

However, these plots suggest that one is learning too slowly from the data. The weights α_n=1/(α+n) decrease to zero too quickly. 
In fact, while, on one hand, one may want small weights α_n for n small, to mitigate the evolution and the consequent sensitivity to permutations, on the other hand, the situation changes after a sufficiently large n. For n large, the underlying random _n will be fairly close to its limit G, and one gets to a situation close to an exchangeable mixture model, with θ_i | G   G. Then, in order to learn from the data about G, the weights α_n should not be too small; that is, they should not converge to zero too quickly. 
The last row of Figure <ref> shows the mixing density estimate g_n 
obtained with weights

    α_n={[     1/(α+n)     n ≤ 500; 1/(α+n)^3/4     n > 500 ].

with α=100. This simple choice is yet very clear in  demonstrating the above intuition. The effect of the ordering is attenuated, and the slower decrease to zero of the weights α_n improves learning. 
Roughly speaking, one uses the first part of the sample to induce some randomness around G_0; in other words, to learn about the prior distribution of G, in a sort of empirical Bayes fashion; then, after a sufficiently large n, the learning process proceeds as in an exchangeable mixture model, with the learned prior law. 

We finally consider inference on the mixing distribution function. 
A slower decay to zero of the α_n gives higher variability of G_n, which is reflected in higher conditional asymptotic variance of its limit G. 
We illustrate such effect in Figure <ref>. The plots show the 95% asymptotic credible intervals for the random mixing distribution function G, given x_1:1000, for data and weights as in panels [2,1] and [3,1] of Figure <ref>. More precisely, we computed the marginal asymptotic credible intervals for G(t_j), for the t_j as evidenced in the plots. The black line is the true mixing distribution used in the simulation. 
For α_n=1/(α+n), with α=100 (left panel), the interval is quite narrow. This is explained by the low variability of G_n, due to the fast decaying weights. With weights α_n as in (<ref>), from  Theorem <ref> we obtain 

    G(t_j) | x_1:n≈(G_n(t_j), V_A,n/2 √(n)) .

The larger credible intervals in the right panel of Figure <ref> reflect the higher asymptotic conditional variance obtained in this case. 



§ PROOFS
 

This section provides the proofs of the previous results. 


Proof of Theorem <ref>. 

(i) For every A∈ℬ(Θ) and every n≥ 0,

    E(G_n+1(A)| X_1:n)=
    (1-α_n+1)G_n(A)+α_n+1 E(P(θ_n+1∈ A| X_1:n+1)|  X_1:n)=G_n(A)

Hence, the sequence (G_n) is a measure valued martingale with respect to the natural filtration of (X_n). 
By Lemma 7.14 in <cit.>, there exists a random probability measure G such that G_n converges a.s. to G, in the topology of weak convergence. 

 (ii) Since, for every A, (G_n(A)) is uniformly bounded, it is a closed martingale. Thus, for every n≥ 0 and k≥ 1,

    P(θ_n+k∈ A| X_1:n)=E(P(θ_n+k∈ A|
    X_1:n+k-1) | X_1:n)=E(G_n+k-1(A)| X_1:n)=
    E(G(A)| X_1:n).

□



Proof of Proposition <ref>. 

Let us denote by _n the sigma field generated by X_1:n, and let  _∞=∨_n _n and Z be a random variable such that Z |_∞∼ G. We have

    E( ∫_Θ h(z) dG(z) | _n.)   =    E( E(h(Z) |_∞) |_n ) = 
    E(h(Z) |_n) 
       =   ∫_Θ h(z) P(Z∈ dz |_n),

where the last equality follows from the fundamental property of regular conditional distributions (see e.g. <cit.>, eq.(2.4)). Noticing that  
P(Z ∈·|_n)= E( P(Z  ∈·|_∞) |_n)= E(G(·) |_n), 
we obtain 

    ∫_Θ h(θ) dG_n(θ) = E(∫_Θ h(θ)dG(θ) |ℱ_n).

Since 

    E(∫_Θ h(θ)dG(θ) |ℱ_n) → 
    E(∫_Θ h(θ)dG(θ) |ℱ_∞)   a.s.

and since G is -measurable, then ∫_Θ h(θ) dG_n(θ)→∫_Θ h(θ) d G(θ) a.s.

To prove the last assertion, notice that, from (<ref>), 
∫_Θ |h(θ)| dG_0(θ)= E(∫_Θ |h(θ)|  dG(θ)). 
Thus, if ∫_Θ |h(θ)| dG_0(θ) < ∞, then the non-negative quantity ∫_Θ |h(θ)| dG(θ) is a.s. finite. 

□



Proof of Theorem <ref>

(i) To prove that (X_n) is c.i.d., it is enough to show that 
P(X_n+2∈ B | X_1:n)=P(X_n+1∈ B | X_1:n), for any n ≥ 0 and any B. This is an immediate consequence of the conditional independence of the X_i, given (θ_n), and of Theorem <ref>. Indeed, denoting by F(·|θ) the distribution corresponding to the density f(·|θ), 
we have 
	
    P(X_n+2∈ B | X_1:n)    =     E(P(X_n+2∈ B |θ_n+2, X_1:n) | X_1:n)
    =∫_Θ F(B |θ) P(θ_n+2∈ dθ| X_1:n) 
       =   ∫_Θ F(B |θ) P(θ_n+1∈ dθ| X_1:n) = P(X_n+1∈ B | X_1:n),

where the third equality follows from (<ref>). 

(ii) Let (t_j, j ∈ J)  be a countable dense set of points in 𝕏. 
By Proposition <ref>, F_G_n(t_j) → F_G(t_j),  
for every ω∈Ω_j with P(Ω_j)=1. Now, let Ω^*=∩_j Ω_j. Being J countable, P(Ω^*)=1,  and for any ω∈Ω^*, F_G_n(t_j) → F_G(t_j) for all t_j. 
For distribution functions, convergence on a countable dense set implies weak convergence. Therefore, we have that, P-a.s., F_G_n converges to F_G in the topology of weak convergence.  Now, notice that F_G is a.s. absolutely continuous, with density f_G. By Theorem 1 in <cit.>, a.s. weak convergence of the predictive measures  to an absolutely continuous random measure  implies that the convergence also holds in total variation. Therefore, P-a.s., F_G_n converges to F_G in total variation, which is equivalent to f_G_nL_1→ f_G. 

(iii) Convergence of the predictive distributions to the random probability measure  F_G implies that (X_n) is asymptotically exchangeable, with directing measure F_G (<cit.>, Lemma 8.2). 

□



Proof of Theorem <ref>

Let ℋ be a a countable, convergence determining class   of bounded continuous functions, and k a positive integer. 
By Theorem <ref>,  (X_n) is asymptotically exchangeable with directing random measure F_G; therefore, for almost all ω=(x_1,x_2,…)

    E(∏_i=1^kh_i(X_n+i)| x_1:n)→∏_i=1^k
     (∫_Θ∫_𝕏h_i(z_i)f(z_i| s_i)dμ(z_i) dG(s_i)(ω)),

for every h_i∈ℋ.
Let ω=(x_1,x_2,…) be fixed in such a way that the above equations hold.   
Then, for every j=1,…,k, the sequence of probability measures (P( θ_n+j∈·| x_1:n)) is tight. It follows that the sequence of joint conditional distributions (P( θ_n+1:n+k∈·| x_1:n)) is tight.
Thus, for every increasing sequence of integers, there exists a subsequence (n_j) and a probability measure Q(ω) such that  
P(( θ_n_j+1,…, θ_n_j+k)∈·|  x_1:n_j)→   Q(ω).  
The proof is complete if we can show that 

    Q(A_1×…×  A_k)(ω)=∏_i=1^k G(A_i)(ω),   ,

because (<ref>) implies that the conditional law of ( θ_n_j+1,…, θ_n_j+k) converges weakly to the product measure G^k. The sequence of random variables 
(∏_i=1^k ∫_𝕏h_i(z_i)f(z_i|θ_n+i)dμ(z_i))
is uniformly bounded and, therefore, it is uniformly integrable. Hence, 

    ∫_Θ^k∏_i=1^k∫_𝕏 h_i(z_i) f(z_i| s_i) dμ(z_i)      P(( θ_n_j+1,…θ_n_j+k)∈(ds_1,…,ds_k)| x_1:n_j)
    
       →∫_Θ^k∏_i=1^k ∫_𝕏 h_i(z_i) f(z_i| s_i) dμ(z_i)  dQ(s_1,…,s_k)(ω),

for every h_1,…,h_k∈ℋ.
On the other hand, by (<ref>),

    ∫_Θ^k∏_i=1^k∫_𝕏h_i(z_i)f(z_i| s_i)dμ(z_i)   P(( θ_n_j+1,…θ_n_j+k)∈(ds_1,…,ds_k)| x_1:n_j)
    
       →∫_Θ^k∏_i=1^k∫_𝕏h_i(z_i)f(z_i| s_i)dμ(z_i) dG(s_1)(ω) … dG(s_k)(ω).

Hence, for every h_1,…,h_k∈ℋ,

    ∫_Θ^k∏_i=1^k∫_𝕏h_i(z_i)f(z_i| s_i)      dμ(z_i) dQ(s_1,…,s_k)(ω)
       =∫_Θ^k∏_i=1^k∫_𝕏 h_i(z_i)f(z_i| s_i)dμ(z_i)dG(s_1)(ω) … dG(s_k)(ω).

Since the model is identifiable,
the class 

    {∫ h(z)f(z|θ )dμ(z) : h∈ℋ}

is separating for 𝒫(Θ). It follows that the class 

    {∏_i=1^k∫ h_i(z)f(z|θ)dμ(z):h_i∈ℋ,i=1,…,k    }

is separating for 𝒫(Θ^k) (<cit.>, Proposition 3.4.6).
Thus, (<ref>) implies (<ref>).

□







Proof of Theorem <ref>. The thesis follows from Lemmas <ref> and <ref>, if we can show that

    sup_nE(∫_Kg_n(θ)^2dλ(θ))<∞  Kλ(K)<∞.

Let K be a fixed compact set, with λ(K)<∞. It holds

    E( ∫_Kg_n(θ)^2dλ(θ)  )   =∫_K E(g_n(θ)^2)dλ(θ)=∫_KE( E(g_n(θ)^2| X_1:n-1))dλ(θ).

By the martingale property of the sequence (g_n), and Jensen inequality, we obtain 

    E(g_n(θ)^2| X_1:n-1)
       = g_n-1(θ)^2 
    E( [1+α_n(f(X_n|θ)/∫_Θ f(X_n|θ') g_n-1(θ') dλ(θ')-1)]^2  | X_1:n-1)
       ≤ g_n-1(θ)^2[ 1+α_n^2(1+∫_𝕏f(x|θ)^2/∫_K f(x|θ')g_n-1(θ')dλ(θ') dμ(x)  )]
       ≤ g_n-1(θ)^2[ 1+α_n^2(1+∫_𝕏∫_K f(x|θ)^2/ f(x|θ')g_n-1(θ')dλ(θ') dμ(x)  )]
       ≤ g_n-1(θ)^2[ 1+α_n^2(1+sup_θ_1,θ_2∈ K∫_𝕏f(x|θ_1)^2/ f(x|θ_2) dμ(x)  )].

Therefore

    E(∫_Kg_n(θ)^2dλ(θ))≤ E(∫_Kg_n-1(θ,)^2dλ(θ))[ 1+α_n^2(1+sup_θ_1,θ_2∈ K∫_𝕏f(x|θ_1)^2/ f(x|θ_2) dμ(x)  )]

Iterating, we obtain

    E(∫_Kg_n(θ)^2dλ(θ))≤∫_Kg_0(θ)^2dλ(θ)∏_i=1^n( 1+α_i^2 M_K)
	
with M_K=(1+sup_θ_1,θ_2∈ K∫_𝕏f(x|θ_1)^2/ f(x|θ_2) dμ(x)  ), which is finite by the assumption (<ref>). By  (<ref>),  
sup_nE(∫_Kg_n(θ)^2dλ(θ))<∞.

□








Proof of Lemma <ref>. 

Since G_n(A)→ G(A) a.s., it remains to show that

    ∫_𝕏 P_G_n(A | x)^2 dF_G_n(x) →∫_{x:f_G(x)≠ 0} P_G(A | x)^2 dF_G(x)   a.s.

By Theorem <ref>, F_G_n converges to F_G in total variation, a.s.; therefore

    ∫_{x:f_G(x) = 0} P_G_n(A | x)^2 dF_G_n(x) 
    ≤ F_G_n({x:f_G(x)=0}) → F_G({x:f_G(x)=0})=0     a.s.

Then, denoting 𝕏_0={x:f_G(x)≠ 0}, we have

    |∫_𝕏_0 P_G_n(A | x)^2 dF_G_n(x) - ∫_𝕏_0 P_G(A | x)^2  dF_G(x)|
    	
       ≤   ∫_𝕏_0 |f_G_n(x) - f_G(x)| dμ(x) + |∫_𝕏_0 P_G_n(A | x)^2(x) dF_G(x) - ∫_𝕏_0 P_G(A | x)^2 dF_G(x)|

The first term converges to zero since f_G_n converges to f_G in L_1, by Theorem <ref>; the second term converges to zero by dominated convergence theorem. Thus, the thesis follows.

□

 


Proof of Theorem <ref>.

For every n≥ 1, let 

    M_n,j={[ √(r_n)(G_n(A)-G_n+j-1(A))                      j≥ 1;                         0                       j=0 ].

and let 

    ℱ_n,j={[ ℱ_n+j-1    j≥ 1;     ℱ_n    j=0, ].

where (_n) is the natural filtration of (X_n). For every n≥ 1, (M_n,j)_j≥ 0 is a zero-mean martingale, with respect to the filtration (ℱ_n,j)_j≥ 0 and ℱ_n,1=ℱ_n⊂ℱ_n+1=ℱ_n+1,1.
Let

    Z_n,j:=M_n,j-M_n,j-1j≥ 1,   U_n:=∑_j≥ 1Z_n,j^2,   Z_n^*:=sup_j≥ 1|Z_n,j|.

The thesis follows from Theorem A.1 in <cit.> if we can show that (Z_n^*) is dominated in L^1 and that (U_n) converges a.s. to V_A.




By definition,
Z_n,1=0 and, for j≥ 2,

    Z_n,j   =√(r_n)(G_n+j-2(A)-G_n+j-1(A))
       =√(r_n)α_n+j-1(G_n+j-2(A) -∫_Af(X_n+j-1|θ)dG_n+j-2(θ)/∫_Θ f(X_n+j-1|θ)dG_n+j-2(θ)).

Since √(r_n)sup_k≥ nα_k→ 0, then  (Z_n^*) is dominated in L^1.

To prove that (U_n) converges a.s. to V_A, we employ Lemma A.1 in <cit.>. To be consistent with the notation therein, let us set b_1=r_1 and, for k≥ 1,

    b_k+1=r_k     a_k=1/b_k^2α_k^2.

Then, we can write

    U_n   =b_n+1∑_j≥ 2α_n+j-1^2(∫_Af(X_n+j-1|θ)dG_n+j-2(θ)/∫_Θ f(X_n+j-1|θ)dG_n+j-2(θ) -G_n+j-2(A)  )^2
       =b_n+1∑_k≥ n+1Y_k/a_kb_k^2,

where 

    Y_k=(∫_Af(X_k|θ)dG_k-1(θ)/∫_Θ f(X_k|θ)dG_k-1(θ) -G_k-1(A)  )^2.
 
It can be proved, proceeding as in Lemma <ref>, that

    E(Y_k|ℱ_k-1)=E( [∫_Af(X_k|θ)dG_k-1(θ)/∫_Θ f(X_k|θ)dG_k-1(θ) -G_k-1(A)]^2|ℱ_k-1)→ V_A   a.s.

as k→∞. 
Moreover,

    ∑_k≥ 1E(Y_k^2)/a_k^2b_k^2<∞,

as ∑_k≥ 1 (a_kb_k)^-2=∑_k≥ 1α_k^4b_k^2<∞ and |Y_k|≤ 1.
Since, by assumption, 
b_n+1∑_k≥ n+1(a_kb_k^2)^-1=r_n∑_k>nα_k^2→ 1,
then, by Lemma A.1 in <cit.>, 
U_n→ V_A a.s. as n→∞. 


□



Proof of Theorem <ref>.
  
We first prove that, for every A∈ℬ(Θ), the conditional distribution of (√(r_n)(G_n(A)-G(A)),V_A,n), given X_1:n, converges to (0, V_A) ×δ_V_A a.s. on the set {ω : V_A(ω) >0}, ω = (x_1, x_2, …). 
To show this,  we use Lemma <ref> and compute the joint characteristic function

    E(exp(i s_1√(r_n)(G(A)-G_n(A))+i s_2V_A,n)| 
    X_1:n)   =E(exp(is_1√(r_n)(G(A)-G_n(A)))| X_1:n)exp(is_2V_A,n)
       →exp(-s_1^2 V_A/2)exp(is_2V_A)

Let now 𝒟 be a countable convergence-determining class of bounded continuous functions for the probability measures on ℝ 
and let 

    D_n=√(r_n)(G(A)-G_n(A)),   
    W_n=1/√(V_A,n)1_(V_A,n>0),
       W=1/√(V_A)1_(V_A>0).
 
Then W_n(ω)→ W(ω) for almost all ω such that V_A(ω)>0. 
By Theorem <ref>, for every h∈𝒟,  

    E(h(D_n)| X_1:n) 1_(V_A>0)→∫ h(x) ϕ(x | 0,1/W^2)dx    a.s.

where ϕ(x|μ, σ^2) denotes the (μ,σ^2) density computed at x.
Since W_n is  a function of X_1:n,  then for every h_1,h_2∈𝒟, 

    E(h_1(D_n)h_2(W_n)| X_1:n) 1_(V_A>0)   = E(h_1(D_n)| X_1:n) h_2(W_n) 1_(V_A>0)
       →∫ h_1(x) ϕ(x | 0,1/W^2) dx h_2(W) 1_(V_A>0)
       =∫ h_1(x_1) h_2(x_2) d((0,1/W^2) ×δ_W)(x_1, x_2)   1_(V_A>0)

Being the class {h_1 h_2:h_1,h_2∈𝒟} a convergence determining class for the probability measures on ℝ^2, then, for every bounded continuous function h,  

    E(h(√(r_n)G(A)-G_n(A)/√( V_A,n))| X_1:n) =
    E(h(D_n W_n)| X_1:n) 
       →∫ h(x W) ϕ(x | 0,1/W^2) dx 
    =∫ h(y) ϕ(y | 0,1) dy .

a.s. on the set {ω : V_A(ω)>0}.
□



Proof of Theorem <ref>
 
Let c_1,…,c_k  be arbitrary real numbers. The sequence  
(∑_i=1^k c_i G_n(A_i))_n≥ 1, is a bounded martingale,  converging a.s. to
∑_i=1^k c_iG(A_i).
Following the same steps as in Theorem <ref>, with ∑_i=1^k c_i G_n(A_i) in the place of G_n(A) and ∑_i=1^k c_iG(A_i) in the place of G(A),
 we obtain 

    P( √(r_n) (∑_i=1^k c_iG(A_i)-∑_i=1^k c_iG_n(A_i)) ≤ t | x_1:n) →Φ(t ; 0, U),   ,

where
U is the a.s. limit of

    U_n:=r_n ∑_j≥ n+1α_j ^2Y_j

with

    Y_j   =[∑_i=1^k c_i(
    P(θ_j∈ A_i | X_1:j)-P(θ_j∈ A_i| X_1:j-1)	
    )  ]^2.

Applying Lemma A.1 in <cit.>, as in the proof of Theorem <ref>, and noticing that

    E(Y_j| X_1:j-1)   =∑_i,i' c_i c_i'[∫_𝕏 P_G_j-1(A_i | x) P_G_j-1(A_i'| x)dF_G_j-1(x)-G_j-1(A_i)G_j-1(A_i')]
       →∑_i,i' c_ic_i'[∫_𝕏 P_G(A_i | x)P_G(A_i'| x)dF_G(x)-G(A_i)G(A_i')]
       =∑_i,i'=1^jc_ic_i'C_A_i,A_i',

we obtain U=∑_i,i'=1^k c_ic_i' C_A_i,A_i'. 
Thus, for every c_1,…,c_k, 

    P( √(r_n) (∑_i=1^k c_iG(A_i)-∑_i=1^k c_i G_n(A_i)) ≤ t | X_1:n) →Φ(t  |  0,∑_i,i'=1^k c_i c_i' C_A_i,A_i' ),
 
for every t. The thesis follows from Cramér-Wold theorem.

□



 Proof of Theorem <ref>.

 Let (_n) denote the natural filtration of (X_n). Consider the _n-measurable spectral decomposition

    C_n(A_1,…,A_k)/r_n=Q_nΛ_n Q_n^T

where Q_n in a k× k orthogonal matrix and Λ_n=diag(λ_1^(n),…,λ_k^(n)). 
Let 
    Y^(n)=Q_n^T[
    [ G(A_1)-G_n(A_1);               …; G(A_k)-G_n(A_k) ]]

and

    Z_i^(n)=Y_i^(n)/√(λ_i^(n))1_(λ_i^(n)≠ 0)+Z̃_i 1_(λ_i^(n)=0),

where Z̃_1^(n),…,Z̃_k^(n) are i.i.d. random variables, independent of ℱ_n, and with (0,1) distribution.  
Then 

    [
    [ G(A_1)-G_n(A_1);               …; G(A_k)-G_n(A_k) ]]=Q_nΛ_n^1/2Z^(n)= C_n(A_1,…,A_j)^1/2/√(r)_n Z_*^(n)

where Z_*^(n)=Q_nZ^(n)≈(0, I), given ℱ_n. 

□



Proof of Proposition <ref> .

With the same notation as in the proof of Theorem <ref>, we can write

    lim inf_n P((G(A_1),…,G(A_k)∈ E_n^(ϵ)|ℱ_n)
       ≥lim inf_n
    P(Z_*^(n)T C_n(A_1,…,A_k)^1/2/√(r_n)( C_n(A_1,…,A_k)+ϵ I/r_n)^-1C_n(A_1,…,A_k)^1/2/√(r_n)Z_*^(n)≤χ^2_1-γ|ℱ_n)
    
       ≥lim inf_n
    P(Z_*^(n)TZ_*^(n)≤χ^2_1-γ|ℱ_n)=1-γ.

□







§ DISCUSSION


For its simplicity and good practical performance, Newton's algorithm is widely used in problems involving hidden variables. We have proposed a novel approach to study its properties, and make users more aware 
of the modeling assumptions implicitly made. We believe that our approach can be useful in other settings, too. Further explicit results on the probability law of the limit mixing distribution G, although difficult, could lead to novel priors for Bayesian nonparametrics on the space of  absolutely continuous  distributions. Modifications of the algorithm could be envisaged (as suggested in the simulation study of Section  <ref>, or by initializing the procedure with exact computations from the DP mixture model)  to control the prior distribution on G.  


A limitation of Newton's algorithm is that it requires to evaluate an integral at each step. An interesting class of recursive algorithms, that allow to overcome this difficulty, has been recently proposed by <cit.>. Multivariate extensions are in <cit.>. Directions of research are extensions of  our study to this class of algorithms, as well as developments for multivariate mixtures and dependent mixture models, possibly exploiting theoretical results on partially c.i.d. sequences (<cit.>).




Acknowledgments This work benefited from discussions with Lorenzo Cappello and Stephen Walker, to whom we express our thanks. The authors acknowledge financial support by the PRIN grant 2015SNS29B.




plainnat   






