{'fraction_non_alphanumeric': 0.03878165397556225, 'fraction_numerical': 0.01941296263502745, 'mean_word_length': 4.328104729331289, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamic Mixture of Experts Models for Online Prediction

A mixture of experts models the conditional density of a response variable using a mixture of regression models with covariate-dependent mixture weights. We extend the finite mixture of experts model by allowing the parameters in both the mixture components and the weights to evolve in time by following random walk processes. Inference for time-varying parameters in richly parameterized mixture of experts models is challenging. We propose a sequential Monte Carlo algorithm for online inference and based on a tailored proposal distribution built on ideas from linear Bayes methods and the EM algorithm. The method gives a unified treatment for mixtures with time-varying parameters, including the special case of static parameters. We assess the properties of the method on simulated data and on industrial data where the aim is to predict software faults in a continuously upgraded large-scale software project.


Introduction
A mixture of experts (ME) model (Jordan and Jacobs, 1994) provides a flexible framework for expressing the distribution of a response variable conditional on a set of covariates. It models continuous or discrete response variables using a finite mixture model with covariate-dependent component models and mixture weights. The component models are commonly referred to as the experts and the mixture weights models as gates. Generally, the experts are regression models with the response's conditional density from the exponential family; see Gormley and Frühwirth-Schnatter (2018) for a concise introduction to ME models and Yuksel et al. (2012, Section IX) for an extensive list of application areas. ME models have been extended in many ways. Hunter and Young (2012) relax the parametric assumption of the component models and propose a semi-parametric inference methodology for mixtures of linear regression models. Villani et al. (2012) extend the component models to density functions outside the exponential family and use Bayesian variable selection in all parts of the model. Jacobs et al. (1997) propose a Bayesian hierarchical model in which the conditional density of the response is expressed as a mixture of ME components. Rasmussen and Ghahramani (2002) propose an infinite mixture of Gaussian process experts. Zeevi et al. (1996) and Carvalho and Tanner (2005a,b) extend the ME model framework to autoregressive time series data where the covariates may include lagged values of the response and they study the properties of the maximum likelihood estimator. Wood et al. (2002) propose a Bayesian inference methodology for ME models and allow smoothing spline components. In longitudinal data applications, Muthén and Shedden (1999) consider fixed and random effect covariates in both components and mixture weights. Quiroz and Villani (2013) propose a ME model with static parameters for longitudinal data but where the subjects are allowed to dynamically change mixture components over time.
While many existing ME models discussed in the previous paragraph are very flexible, they are still too restrictive for many applications, in particular when: i) the data is an irregularly spaced time series or arrives in batches; ii) the conditional density of the response tends to change over time. Our application to predicting software faults, in Section 4, is one such typical example taken from industry. Software is generally upgraded at irregular times -a week may pass without any release followed by two or three releases in the next week -depending on the bugs reported in previous versions, the complexity of the added features in the release, or business related factors. Also, the distribution of the response variable naturally changes as the software matures: the developers in the project change over time, user behavior changes, and new technologies emerge which demand adaptation. All these issues make standard ME models for time series data impractical.
We extend the class of ME models to dynamic mixture of experts (DME) models to propose an online (real-time) predictive model with time-varying parameters. The new class is particularly suited for unstructured streaming data, but may also be used for equally spaced time series with time-varying distributions. The proposed DME models have the general form of varyingregression coefficient models (Hastie and Tibshirani, 1993); the regression coefficients in both the mixture components and the mixture weights are allowed to vary over time through a latent process. The latent process can be modelled as a discrete first order Markov process with static (Liehr et al., 1999;Kohlmorgen et al., 2000) or time-varying (Wang et al., 2003) transition probabilities. Alternatively, it can be modelled through random walk processes (West et al., 1985;Fahrmeir et al., 2004;Fahrmeir and Kneib, 2011) as in Section 2.1. The density for the response variable is modelled conditional on the value of the latent state. This allows the model to adapt locally to any abrupt changes through time (Fahrmeir et al., 2004). Also, it facilitates designing an adaptive online (realtime) Bayesian predictive inference in which the uncertainty is updated sequentially as new data arrive and predictions are computed based on the recently updated posterior density -the online posterior; see Section 3.
Inference in mixture models is challenging. One challenge is the identifiability issue caused by the invariance of the likelihood under permutation of component labels and parameters. Constraints on the regression coefficients have been imposed to identify ME models, see (Jiang and Tanner, 1999) and Appendix B. However, identifiability is of less importance in applications where prediction is the main objective as the predictive distribution is unaffected by label switching (Geweke, 2007). Another important issue in mixture models is the number of components. One strand of the literature use nonparametric Bayes approaches with infinite mixtures typically modeled by stick-breaking processes (Hjort et al., 2010), or reversible jump MCMC (Richardson and Green, 1997) to obtain the posterior over the number of mixture components. Another line of research uses model comparison methods to select the optimal number of components (Celeux et al., 2019;Geweke and Keane, 2007). In addition to these issues, the structure of DME models brings additional challenges. DME are often richly parameterized with time-varying regression coefficients in both the components and the mixture weights. This leads to a high dimensional and complex target posterior that can only be properly explored by carefully designed numerical methods.
We propose online inference based on Sequential Monte Carlo (SMC) methods (Del Moral et al., 2006;Doucet et al., 2000) to address these issues. The performance of SMC methods depends on the proposal/importance distribution from which parameters are sampled. Given the complexity of DME models, off-the-shelf SMC methods based on commonly used proposal distributions such as the bootstrap filter will perform poorly. Our main contribution is a proposal distribution tailored to the class of DME models. The potentially high-dimensional regression coefficients influence the conditional density only through the low-dimensional linear predictors transformed by the link functions. This makes it possible to combine the linear Bayes method (West et al., 1985) and ideas from the expectation and maximization (EM) algorithm principle (Dempster et al., 1977) to build the proposal distribution efficiently. Our algorithm builds on the marginal particle filter algorithm (Klaas et al., 2005) for online (real-time) prediction; offline inference using SMC algorithms can be done following ideas in Munezero (2021), where the proposed methodology can be considered as an integrated part of the particle smoother.
Standard inference for mixture models requires posterior draws of component allocation indicators or, in the approach of Carvalho et al. (2010), to keep track of the number of allocations for each component and conditional sufficient statistics for the mixture component parameters. Our methodology expresses the likelihood in its marginal form, which does not require sampling component indicators, hence reducing considerably the dimension of the target posterior. The posterior's dimension is reduced even further by using the discount factor approach (West et al., 1985) that recursively estimates the innovation variance of the states at a particular time point as a function of a discount factor 0 < α < 1 and the filtered information from the most recent posterior. This means that the inference requires keeping track of only the regression coefficients through time.
The smoothness of the regression coefficients evolution is controlled via α which allows building static and dynamic models in a unified way just by changing the value of α. Liu and West (2001) suggest that models with 0.95 ≤ α < 1 are essentially static, and those with α < 0.95 are dynamic, although this rule of thumb may vary depending on applications. Inference on the discount factor and the number of components is performed using the log predictive score, a marginal likelihood-based model comparison criteria (Geweke and Keane, 2007;Villani et al., 2009).
The proposed methodology allows writing general computer code where a user can easily add a new model by supplying the first and second derivative of the component densities with respect to the linear predictors, which are the only key arguments of our procedure. Recent advances in automatic differentiation (Baydin et al., 2018) even removes the requirement of computing derivatives analytically.
The rest of the paper is organized as follows. Section 2 introduces the dynamic mixture of experts model and the prior process. Section 3 presents the SMC algorithm based on a proposal distribution from linear Bayes theory. Section 4 presents an industrial application to online prediction of faults in a large-scale software project, where allowing parameters to evolve over time considerably improves predictive performance. Section 5 explores the properties of the inference method on simulated data. The final section concludes.

Dynamic mixture of experts
Let D j = (y j ,x j ) represent data from a time dependent process observed at different time points j = 1, . . . , J, where y j denotes the univariate response variable andx j = (x (1) j , . . . ,x (P) j ) is a Pdimensional covariate vector. The D j may contain only one observation as in standard time series applications, or it may be a data batch containing several observations as in the software upgrade process described in Section 4. We propose the dynamic mixture of experts model for online (real time) prediction of y j given the value of the covariatex j ; z j and x j are subsets of x j of dimensions Q and P respectively. The λ jk (x j ) and ω jk (z j ), k = 1, . . . , K, are time-varying covariate-dependent parameter and mixture weight functions of the k th expert model respectively, λ j = (λ j1 (x j ), . . . , λ jK (x j )) and ω j = (ω j2 (z j ), . . . , ω jK (z j )). The covariates in the mixture weights can be distinct from the covariates in the experts.
Here, the experts represent the component models in the mixture (2.1) and they depend on the structure of the response variable; they are typically density functions from the exponential family, e.g. Gaussian if y j is continuous, or Poisson, binomial or negative binomial for count data, or multinomial if y j is categorical. However, as in Villani et al. (2012), we allow the component models to be any well-behaved density functions, not necessarily limited to the exponential family, and the model parameter may be multidimensional with each of its components connected to the covariates through its own link function. By well-behaved densities we mean densities that are twice differentiable with respect to the parameters, and that satisfies the nondegeneracy condition in Jiang and Tanner (1999, Condition 1) so that the mixture of experts model is identified. Jiang and Tanner (1999) show that this condition is fulfilled for Poisson components and Section B in the Appendix give a similar result for generalized Poisson components, which do not belong to the exponential family and are used in the empirical application in Section 4.
The component model parameters λ jk = λ jk (x j ), k = 1, . . . , K are connected to their linear predictors through a link function g as jk ) . For component models with more than one parameter, Eq. 2.1 can be extended by linking each parameter to its own linear predictor; see Villani et al. (2012). Furthermore, the mixture weights depend on the covariate z j , through the multinomial logit link function jk ) . Following standard practice we set ψ j1 = 0 for all j in (2.3) to identify the model (Jiang and Tanner, 1999, Remark 1). In the following, we refer to β jk and θ jk as the regression coefficients in the component distributions and mixture weights, respectively, and to λ j = (λ j,1 , . . . , λ j,K ) and ω j = (ω j2 , . . . , ω jK ) as the model parameters.
Mixture models are well-known to suffer from label switching, i.e. invariance under permutations of the components. A common approach is to impose order restrictions on parameters of the mixture components Jiang and Tanner (1999), either before running MCMC or SMC to sample from the posterior, or by re-ordering the posterior draws after the sampling (Stephens, 2000). Alternatively, in the Bayesian framework, the identifiability problem is addressed by designing informative priors (Malsiner-Walli et al., 2017). However, our interest here is on predictive inference, where label switching is not a concern (Geweke, 2007).

Prior process
To simplify notation, we stack all the regression coefficients at time j into one vector γ j = (β j , θ j ) , where β j = (β j1 , . . . , β jK ) and θ j = (θ j2 , . . . , θ jK ), and the linear predictors for all components into ρ j = (η j , ψ j ) . The prior for the γ j is a random walk with a predefined initial distribution p(γ 1 ), which allows it to vary over time. This prior process is commonly applied in dynamic models as a way of penalizing the regression coefficients from high fluctuations and avoiding overfitting (Fahrmeir and Kneib, 2011;Fahrmeir et al., 2004). In some applications it is sufficient to set U j = U, which is a special case of (2.5). However, it is more useful for online inference to let U j change over time as it allows to update the prior with historic data recursively as more data batches are observed. Fully Bayesian inference requires a prior for each U j . Common priors are: i) an inverse-Wishart density for a full matrix U j (Gamerman, 1998), ii) an inverse-gamma density (Fahrmeir et al., 2004) or a random walk process (Lang et al., 2002) on the elements of a diagonal U j . An alternative to placing a prior on each U j is to approximate each U j recursively using the discount factor approach in West et al. (1985). Let C j denote the posterior covariance of γ j and set U j = (α −1 − 1)C j−1 for a given discount factor 0 < α < 1. A value of α close to one shrinks U j towards zero, leading to very little variation in γ j over time; a value of α close to zero gives the regression parameters more flexibility and allows the model to adapt well to local fluctuations in the parameter; for instance, change points or level shifts in the parameter.
The discount factor approach has some advantages compared to a fully Bayesian approach. It is computationally much quicker as it avoids extra simulations from the posterior of U j . The discount factor conveniently controls the smoothness of the parameter evolution through time with a single parameter, and it allows building static and dynamic models in a unified way just by changing the value of α. Following Liu and West (2001), models with .95 ≤ α < 1 are essentially static, and those with α < .95 are dynamic. We use this approach in Sections 4 and 5.
Our inference methodology applies also to the case of a fully Bayesian approach where U j is estimated in an additional step using particle Markov chain Monte Carlo (Andrieu et al., 2010) and SMC2 (Chopin et al., 2013) methods which allow inference in models with both fixed and time-varying (latent) parameters. It can also be used in the online parameter learning methodology of Carvalho et al. (2010).

Inference, prediction and model comparison
The state space model in Section 2 enables us to exploit the vast literature (Gordon et al., 1993;Pitt and Shephard, 1999;Doucet et al., 2000Doucet et al., , 2006Doucet and Johansen, 2009;Klaas et al., 2005) available on sequential Monte Carlo (SMC). SMC methods are particularly appropriate for sampling from the online posterior and real-time predictive distributions. The present model often has many parameters, and off-the-shelf SMC algorithms with simple proposal distributions like the bootstrap filter (Gordon et al., 1993) will therefore perform poorly. This section describes our proposed algorithm for sampling from the online posterior using a particle filter tailored specifically to the class of dynamic mixture of experts models. We also present the model comparison criteria used to select the number of mixture components and the discount factor.

The marginal particle filter approximation of the online posterior distribution
The target density is the online posterior p γ j | D 1:j updated sequentially in time using a prediction step followed by a measurement update step using Bayes' theorem to make prior-to-posterior updates. The function f j (·) is the response density defined in (2.1), D 1:j denotes the data observed until time j, and p γ j | D 1:j−1 is the prior updated with all historic data observed before the data batch D j . Note that contrary to (2.1), now f j (·) is parametrized in terms of only the regression coefficients γ j ; this is because all quantities required in the inference methodology discussed later are expressed in terms of γ j only.
We are interested in the online predictive distribution p y j |x j , y 1:j−1 which only depends on the filtering density up to time j − 1 (Doucet et al., 2000). However, the challenging part of the sequential inference in (3.1) -(3.2) is that the integral in (3.1) is only tractable for linear Gaussian models (West et al., 1985;Gordon et al., 1993). One way to sample from (3.2) is to use a particle filter algorithm. The particle filter is very attractive for real-time predictions; it allows to sample from intractable distributions and it does not require a scan of the full dataset every time a new observations becomes available. We use the marginal particle filter of Klaas et al. (2005) to generate a set of particles {γ m j } M m=1 associated with the importance weights {w m j } M m=1 . Given the particle sample, any posterior expectation The estimator (3.4) converges to E h(γ j ) as M → ∞ under some weak assumptions stated in Geweke (1989); see Doucet et al. (2001), and Chopin (2004) for more results on the convergence of the particle filter. The marginal particle filter proposes particles from the proposal distribution q(γ j | D 1:j ) and computes the importance weights as The estimation of the density in (3.1) follows Eq.(3.4); i.e, Notice that the importance weights in (3.5) depend on the likelihood f j (y j |x j , γ j ) expressed using the covariates and the regression coefficients rather than the model parameters as in (2.1); the inference requires keeping track of the regression coefficients only.
Clearly the proposal density plays an important role. A proposal density which is inconsistent with the target posterior may lead to the particle degeneracy: The importance weights of only a few particles tend to be substantially different from zero, leading to very few effective samples. To mitigate this degeneracy issue, particles with low weights are discarded and replaced by copies of the particles with high weights. Various strategies for resampling particles are available in the literature (Gordon et al., 1993;Liu and Chen, 1998;Carpenter et al., 1999;Fearnhead and Clifford, 2003) and Douc and Cappé (2005) compare some of these resampling schemes.
In the next section we use linear Bayes methods (West et al., 1985) to construct a proposal q(γ j | D 1:j ) that is tailored to the true posterior, which is crucial for particle methods in highdimensional parameter spaces. West et al. (1985) develop a linear Bayes method (Goldstein and Wooff, 2007) for dynamic generalized linear models with recursions for the posterior mean and covariance over time, making no assumptions on the distributional form of the posterior. Ravines et al. (2007) use these recursive moments to design a multi-move proposal for MCMC targeting the joint smoothing posterior in dynamic generalized linear models. We combine the linear Bayes method in West et al. (1985) with ideas from the EM algorithm Bishop (2006, ch.9) to design a proposal distribution q γ j |D 1:j targeting the filtering density p(γ j | D 1:j ) in dynamic mixture of experts models. The proposed method allows general mixture components outside the exponential family with any twice differentiable link function. Similar to Eq. (2.8) in West et al. (1985), we can write the joint posterior of the regression coefficients and the linear predictors as

A computationally fast proposal distribution for high-dimensional marginal particle filters
where we recall that ρ j = (η j , ψ j ) contains the linear predictors in all components and mixture weights. The second factor in (3.7) does not condition on D j since γ j only enters the likelihood function through the scalar-valued linear predictors in each component, η jk = x j β jk and ψ jk = z j θ jk for k = 1, . . . , K. Our proposal is tailored to the posterior p γ j |D 1:j by using the following steps: 1. Approximate the prior p(γ j |D 1:j−1 ) using a Gaussian with mean and covariance computed from particles at time j − 1.
3. Propose from q(γ j |D 1:j ) = N(µ j , H j ), where µ j and H j are obtained from the law of iterated expectation and law of total variance on (3.7) using a Gaussian approximation of p(ρ j |D 1:j ).
To give the details of the three steps, define η j := X j β j , where X j = I K ⊗ x j and ψ j := Z j θ j , where Z j = I K ⊗ z j ; hence, we can compactly write ρ j = W j γ j where γ j = (β j , θ j ) , ρ j = (η j , ψ j ) and We can use particles from time step j − 1 to approximate γ j |D 1:j−1 ∼ N(γ j , Σ γ j ), wherē and then obtain the mean and covariance of the second factor of (3.7) by conditioning this distribution on the linear constraints ρ j = W j γ j (Rue and Held, 2005, eq. 2.28-2.29) yielding

the proposal is
q γ j |D 1:j = N µ j , H j with moments obtained from applying the law of iterated expectations and the law of total variance to (3.7), It remains to compute E ρ j (ρ j |D 1:j ) and V ρ j (ρ j |D 1:j ). A second order Taylor expansion of log p(ρ j |D 1:j ) aroundρ j leads to the following approximations (Doucet et al., 2000): (3.11) Letting π jk = log ω jk f jk y j |λ jk , the gradient can be computed by direct calculation where Pr(s j = k|D 1:j ) ∝ ω jk f jk (y j |D 1:j−1 , λ jk ) are the posterior probabilities of the observation y j coming from component k (see Bishop, 2006, ch. 9.3 for similar expressions for the EM algorithm). Similarly, the Hessian is, Pr(s j = k|D 1:j )∇∇ ρ j π jk − Σ −1 ρ j Note that the component parameters η jk and ψ jk enter additively in log ω jk f jk y j |λ jk ; therefore, their gradients can be computed separately.
If the batches D j contain several observations, then µ j and H j can be computed by iterating the procedure described above over the observations in the batch; see Gamerman (1991) for a similar approach. Starting with the first observation, we proceed through the following iterations: 1. Compute µ (i) j and H (i) j from Eq. 3.9 and Eq. 3.10.

Setγ
3. Return to step 1 until the last observation in the batch.

Model comparison and prediction
Our model depends on the choice of the number of mixture components K and the discount factor α. We propose to infer those quantities using a sequential version of the marginal likelihood (Doucet et al., 2000) p (y 1:J ) = p (y 1 ) J ∏ j=2 p y j |y 1:j−1 , (3.12) where p y j |y 1:j−1 = f j y j |x j , γ j p γ j | D 1:j−1 dγ j . (3.13) Given a sample of M particles {γ m j−1 } M m=1 and the corresponding importance weights {w n j−1 } N n=1 , the predictive distribution (3.13) is approximated aŝ where γ m j are generated from the transition distribution p(γ j |γ m j−1 ). Different predictive scores are defined as functions of (3.12). One particular example is the log predictive score where J * ∈ [1, . . . , J). The LPS is generally sensitive to the initial distribution of the parameters (Villani et al., 2009). We therefore use the last J /2 data batches to compute the LPS for the models in Section 4; i.e J = J /2. We assume that the particle approximation to the marginal likelihood should be stable after j = J /2. Computing the LPS for different combinations of the number of mixture components K and the discount factor α makes it possible to select good values for these model specification parameters.

Predicting faults in large-scale software projects
Large-scale industrial software projects are continually upgraded to fix bugs and/or to add new features. The upgrades are generally at irregular times: in one week we may observe one release, and in the next, two or three releases depending on several factors such as the amount and severity of the bugs reported in previous versions, the complexity of the new features added to the software, and other business-related factors. Other key factors include the human interaction with the software and the technology evolution. The developers, the user behavior and technologies change over time. The dynamic mixture of experts model (2.1) is appropriate in this case. Different mixture components allow us to model the unknown variations/changes in the human interaction with the software and the time-varying parameters enable the model to adapt to the changes over time.
As the response variable is the number of faults y t reported on the upgrade created at time t, we propose a dynamic mixture of Poisson experts. Here each expert is a Poisson regression model with a covariate vectorx t selected from six code complexity metrics that measure changes made in the source code. The metrics include: i) The number of commits (NC) which represents the number of modifications done from the previous to the current version, ii) the number of changed modules (CM), iii) the number of faults corrected (NFC) per line of code which is the ratio of the total number of faults corrected and the total number of code lines excluding comments, iv) the proportion of files written in C++ (CF), v) the proportion of files written in Java (JF), and vi) the file complexity (FC). The latter is a score calculated based on the number of control flows in the code, e.g. if, for and while loop statements.
The aim is to build an online prediction model for the number of faults in a planned upgrade release. We use a software trouble reports data set from a large-scale project at a major telecom company; the dataset contains a history of 1800 upgrades that were created during a period of 650 days (roughly 21 months). All covariates, excluding the CF, JF and NFC are integers ranging from zero to a value up to six order of magnitude. Therefore, to reduce the scale variations, we apply the log(1 +x t ) transformation to the integer complexity metrics; after this transformation the highest value is no greater than 15.
To make it tractable to deal with the irregular times of fault reports, we partition time into short contiguous intervals [τ 0 , . The partition of time induces a partition of the original data into a sequence of batches D j = y j ,X j which collect data for all upgrade packages created within the time interval t ∈ [τ j−1 , τ j ), j = 1 . . . , J. Batch D j contains N j data points, where y j = (y 1j , . . . , y N j ,j ) is a vector of the response observations in the batch andX j = (x 1j , . . . ,x N j ,j ) is a vector of covariatesx ij for the data point i = 1, ..., N j .  The time is partitioned into 30 days-long intervals, which leads to 21 intervals in total. Experimentation with intervals lengths of one week, two weeks and three months did not improve the LPS. We also assume the initial distribution γ 1 ∼ N(0, I), where I is the unit diagonal matrix; see Appendix A.2 for details. Table 1 compares different fitted models based on their LPS. The table displays various dynamic models, with discount factor α = 0.5, and their static versions, where α = 0.99. The models in the table have different variables in the component models and the number of commits (NC) as the only covariate z in the mixture weights. To select z, we fix the covariates in the component models toX and, starting from z =X, we eliminate variables in z systematically based on the LPS.  Table 1 shows that dynamic models outperform static models, with a difference in LPS more that 40 for single component models and well above 100 for several of the multicomponent models. Also, there is a very large jump in LPS when going from one to two components, in particular for the dynamic versions. While two components seem to be sufficient for the dynamic models, the static models require more components and covariates. The dynamic model CM+FC with two components seems to perform well in terms of LPS since adding more complexity gives no significant increase in LPS. This model is therefore selected for further analysis.
To illustrate that our algorithm can be used also for models outside the exponential family we fit a single component dynamic generalized Poisson model (Famoye and Singh, 2006) using CM and FC as covariates in both the mean and dispersion functions; see Appendix A.3. Figure  4.2 displays the predictive distribution for the one and two-component versions of the selected CM+FC dynamic model, and the dynamic generalized Poisson model at three time points: j = 2, j = 10 and j = 21. The predictive distribution at the time point j is constructed using the posterior at the previous interval j − 1 and the batch D j as test set. To investigate the efficiency of the proposed SMC inference methodology, we fit the selected two-component CM+FC dynamic model using a particle filter with 1000 particles. This is at least an order of magnitude smaller than what can easily be afforded in real applications, but is used here to investigate how much the inferred predictive distribution varies over 100 independent runs with different seeds. Figure 4.3 shows that this variability is small; the figure also includes the predictive distribution from a single run with 100, 000 particles to represent the ground truth. This shows that the proposed method is very efficient and even a small number of particles gives an adequate numerical precision for most applications. the 95% intervals of the predictions from 100 independent iterations of the particle filter with 1000 particles. The orange curve represents the prediction obtained from a particle filter with 100, 000 particles.

Simulation study
We perform several simulation experiments to study the performance of the proposed inference methodology on data generated from both static and dynamic data generating processes (DGPs).

Simulation experiments
The simulation experiments simulate data from the five data generating processes summarized in Table 2. Models M 1 , M 3 , and M 3 all assume a batch data structure, which is common in industrial applications, where data are observed at irregular time points and aggregated into batches. The only source of time variation is in the parameter evolution; the parameter stays constant within a batch but may change across batches. Models M 4 and M 5 are mixtures of autoregressive experts Tanner, 2005b, 2007) with constant (M 4 ) and time-varying parameters (M 5 ). The latter two models are pure time series models with one data point observed at equidistant times, and where the response depends on its lagged values. We follow Carvalho and Tanner (2007) and use log(y t−1 + 1) as lagged values, and parametrizing the autoregressive parameter as 1 − exp(β), so that β can vary freely while at the same time ensuring that the process is stationary for every parameter value.
For each DGP, 50 datasets of 1000 observations are generated. For M 1 to M 3 , data are generated sequentially over 10 time intervals, having 100 observations within each interval. For M 4 and M 5 a time series of length 1000 is generated. The first half of the data is used for training and the last for validation; all model comparisons are based on LPS values computed on the test set -the last half of generated data. Table 2: All data generating processes used in the simulation study. For M 1 , M 2 and M 3 the covariates x ij and z ij are iid U(−1, 1).

Inference of the number of mixture components and the discount factor
The number of mixture components/experts K and the discount factor α discussed in Section 2 are unknown. Inference of these quantities is a research area. Here, we use LPS to assess the performance of the proposed methodology on the inference of these quantities. Several models with k = 1, 2, 3 Poisson components and α ∈ {0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.99} are fitted to data generated from each of the M 1 , M 2 and M 3 DGPs, and LPS is used to select the best model. The aim is to see if the proposed inference methodology is able to identify the underlying data generating process. Figure 5.1 displays the selection frequency of K and α for all fitted models. For M 1 and M 2 , the most frequently selected model is the single component Poisson model with α = 0.99, and α = 0.4 respectively. While, for M 3 , it is the model with α = 0.6, and not the correct two components mixture model. This slight overestimation is not surprising as LPS is often observed to have a tendency to be generous with the number of components in a mixture without having a large impact on the final predictive density, see e.g. Villani et al. (2012).

Comparing static and dynamic models
The data generating process is generally unknown in real applications and the usual strategy in modeling the data is to fit static models. It is therefore interesting to evaluate how fitting a dynamic model would differ from its static version in the cases where the true data generating process is static or dynamic.
We first consider the M 1 , M 2 and M 3 data generating processes. Figure 5.2 compares the performance of the model i) with K = K opt and α = α opt , where K opt and α opt are the values chosen from LPS and ii) the corresponding static model with K = K opt and α = 0.99. The figure shows boxplots of the difference in the LPS values in the validation set for both models. For M 1 the average LPS difference between the selected and the static models is around zero, which shows that the dynamic model does not overfit on static data. On the other hand, for the two dynamic data generating processes, M 2 and M 3 , the dynamic model selected in the validation step clearly outperforms the static model and the difference in LPS increases with the number of components. Consider now the mixture of autoregressive Poisson experts model, M 4 and M 5 . Static models are fitted using the MCMC algorithm in Villani et al. (2012) with 10000 MCMC iterations and dynamic models are fitted using 1000 particles. We partition the data into batches when running our algorithm, which also allows us to investigate the effect of the chosen batch size. Three different batch sizes are compared: 10, 25 and 50. The MCMC algorithm is also updated sequentially at each batch for comparability and for reducing computing times. The dynamic models are trained with different discount factors α ∈ {.3, .4, .5, .6, .7, .8, .9, .99} and the LPS is used to select the best model in the validation step.  the selected dynamic model and the static model is around zero, which again shows that the dynamic model does not overfit on static data. On the other hand, for the two dynamic data generating processes, the dynamic model clearly outperforms the static model. Also, one can note that the data partition has minimal effect as there is not much variation in the LPS differences for the different batch partitions.

Efficiency of the linear Bayes proposal
The efficiency of particle filter algorithms is generally assessed based on the effective sample size where w m j are the importance weights computed at interval j. To assess the performance of the linear Bayes proposal strategy, we compare it with the local linearisation proposal strategy (Doucet et al., 2000); the state-of-the-art method of constructing proposal densities which approximates the target density (3.2) by a linear Gaussian distribution obtained from a second order Taylor expansion of the target density with respect to the regression coefficients γ j , for j ∈ (1, · · · , J). The main difference between these two methods is that the local linearisation method approximates the target density by a Gaussian density without the intermediate step of updating the linear predictors. Figure 5.4 compares the effective sample size per second generated by the two strategies as a way of comparing their efficiency and computation time. The results are based on data simulated from M 3 and a posterior distribution approximated by 1000 particles and α = 0.5. The figure shows that the linear Bayes proposal generates an effective sample size that is on average 10% higher than the local linearisation. Both methods are quite fast; their computation time on a simple windows laptop with intel core i5 processor is less than 3 CPU minutes.

Conclusions
We introduce a general class of dynamic mixture of experts models for online predictions; the model allows the regression coefficients in each mixture component and weight to vary over time. The component models can be essentially any density function, not necessarily limited to the exponential family.
We propose an efficient SMC algorithm for sequential inference and online prediction that is tailored to handle the proposed model class with potentially high-dimensional parameter spaces. The algorithm handles models with static and dynamic parameters in a unified way.
The model is applied to online prediction of the number of faults in a continuously upgraded large-scale industrial software project. We show that allowing the parameters to evolve over time greatly improves the model's predictive performance. A simulation study documents that the proposed model selection procedure is i) effective in reducing flexibility when data comes from a static single-component model, ii) able to fit data from multi-component models with time-varying parameters, and iii) it is fast and generates an effective sample size rate that is superior to the stateof-the-art particle filter which uses a proposal density designed via the local linearization of the target density. Muthén, B. and K. Shedden (1999). Finite mixture modeling with mixture outcomes using the EM algorithm. Biometrics 55 (2)

A.1 Gradient and Hessian for general dynamic mixture of experts models
This section provides details on the expressions of the gradient and the Hessian of the model (2.1) discussed in Section 3.2. Let π jk = log ω jk f jk y j |λ jk as in Section 3.2, the gradient The Hessian is Pr(s j = h|D 1:j )∇ ρ j π jh ∇ ρ j π jk It can easily be shown that the first term in the expression above is zero; hence ∇∇ ρ j log p(ρ j |D 1:j ) = Pr(s j = k|D 1:j )∇∇ ρ j π jk − Σ −1 ρ j .

A.2 Poisson experts
This appendix provides details on the mixture of experts model with Poisson components fitted to the software reports data. The mixture of experts model with Poisson components for the batch D j , j = 1, . . . , 21, has the form To compute the gradient and Hessian required in the proposal density, we define (omitting the index i) The first derivatives are ∂l k ∂η jk = y j − exp η jk , ∂l k ∂ψ jk = 1 − exp ψ jk 1 + ∑ K k=2 exp ψ jk and the second derivatives are

A.3 Generalized Poisson experts
This section delineates the one-component dynamic generalized Poisson model fitted in Section 4.
To simplify the notation we omit the batch index j. The generalized Poisson model is of the form (Famoye and Singh, 2006) f (y|µ, ϕ) = µ 1 + ϕµ where µ > 0 is the mean and ϕ > 0 is the overdispersion parameter. The mean and the dispersion parameters are connected to covariates via a log links µ = exp(x β) and ϕ = exp(x θ) respectively.
Let l = log f (y|µ, ϕ), the first derivatives are (1 + ϕy) (1 + ϕµ) 2 , and the second derivatives are (1 + ϕy) 2 (1 + ϕµ) 3 B Identifiability Jiang and Tanner (1999) prove that mixtures of generalized linear model experts are identifiable if four conditions hold: i) the experts are irreducible, i.e. no pair of experts have identical parameters; ii) there is an ordering of the parameters to avoid so-called label switching of components; iii) the parameters in mixing functions are set to zero for one of the components; and iv) a certain nondegeneracy condition holds that precludes exact linear combinations of the expert densities. The first condition can be assumed to always hold, otherwise we can just collapse identical experts to a single one and reduce the number of mixture components (Jiang and Tanner, 1999). Condition ii) is rarely explicitly imposed in mixture models since it complicates inference, and unrestricted inference will therefore returns one of the K! identical modes; this is acknowledged in Jiang and Tanner (1999) who recommend in Remark 1 to report the mode corresponding to the order restricted parameters for interpretation. Condition iii) is explicitly imposed in our models by zero restrictions, as is commonly done in mixture of experts models. Condition iv) is a technical condition that should be checked for each distributional family on a case by case basis, but is instead often silently assumed to hold. Jiang and Tanner (1999) prove that several commonly used distributions satisfy the nondegeneracy condition, including the Poisson. In fact, only for binomial experts with the number trials smaller than 2K − 1 do Jiang and Tanner (1999) find that the nondegeneracy does not hold, and even then they conjecture that it will hold "for almost all parameters". Jiang and Tanner (1999) prove their results for generalized linear regression (GLM) components, i.e. for densities in the one-parameter exponential family with a scalar dispersion parameter. However, their proofs only rely on properties of the translation and permutation groups acting on densities, and do not use specific properties of exponential families. Their results therefore also apply to a mixture of GLM-type experts with a single parameter depending covariates, but with a density that may be outside of the exponential family; an example of such a model is the generalized Poisson regression where the mean depends on covariates via a linear predictor through a link function, and the overdispersion parameter is a constant. We will now show that this model satisfies the nondegeneracy condition and a mixture of such experts is therefore identified. For the more general model used in Villani et al. (2012) and in Section 4, where the overdispersion parameter is also allowed to depend on covariates, the issue of identification is not yet resolved. We conjecture however that this model is also identified based on the discussion in the previous paragraph and on the empirical results in Villani et al. (2012) where no convergence issues were encountered when using that model.
To prove that a mixture of generalized Poisson regression experts satisfies the nondegeneracy condition we follow the same technique as in the proof of Lemma 3(a) in Jiang and Tanner (1999) for the identifiability of mixture of Poisson experts. The nondegeneracy condition (Condition 1 in Jiang and Tanner (1999)) is that { f (y|µ k , ϕ k )} 2K k=1 are 2K linearly independent functions of y for any 2K distinct pairs (µ k , ϕ k ), for k = 1, . . . , 2K, i.e. that We can write (y −1 + ϕ k ) y−1 µ k 1 + ϕ k µ k y exp −µ k (1 + ϕ k y) 1 + ϕ k µ k = exp a k + b k y + (y − 1) log(y −1 + ϕ k ) , (B.3) where a k = −µ k /(1 + ϕ k µ k ) and b k = log(µ k /(1 + ϕ k µ k )) − µ k ϕ k /(1 + ϕ k µ k ). Now, as y → ∞ the expression in (B.3) behaves as exp(ã k +b k y) whereã k = a k − log ϕ k andb k = b k + log ϕ k . Similar to the proof of Lemma 3(a) in Jiang and Tanner (1999) we therefore have that as y → ∞ the decay/explosion rates of exp(ã k +b k y) are different for each k since the µ k and/or ϕ k are different. Hence since the equality in (B.1) needs to hold for all y ∈ {0, 1, . . .}, this can only be true for c 1 = c 2 = . . . = c 2K = 0. The mixture of generalized poisson experts therefore satifies the nondegeneracy condition.