Poisson-Gamma Dynamical Systems with Non-Stationary Transition Dynamics

Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 06:46:41

Abstract: Bayesian methodologies for handling count-valued time series have gained prominence due to their ability to infer interpretable latent structures and to estimate uncertainties, and thus are especially suitable for dealing with noisy and incomplete count data. Among these Bayesian models, Poisson-Gamma Dynamical Systems (PGDSs) are proven to be effective in capturing the evolving dynamics underlying observed count sequences. However, the state-of-the-art PGDS still falls short in capturing the time-varying transition dynamics that are commonly observed in real-world count time series. To mitigate this limitation, a non-stationary PGDS is proposed to allow the underlying transition matrices to evolve over time, and the evolving transition matrices are modeled by the specifically-designed Dirichlet Markov chains. Leveraging Dirichlet-Multinomial-Beta data augmentation techniques, a fully-conjugate and efficient Gibbs sampler is developed to perform posterior simulation. Experiments show that, in comparison with related models, the proposed non-stationary PGDS achieves improved predictive performance due to its capacity to learn non-stationary dependency structure captured by the time-evolving transition matrices.

Main paper: 
 introduction: 
 In recent years, there has been an increasing interest in modeling count time series. For instance, some previous works [1, 2, 3] are concerned with how to learn the evolving topics behind text corpus (frequencies of words) over time. Some works [4, 5, 6, 7] try to predict global immigrant trends underlying international population movements. Count time series are often overdispersed, sparse, high-dimensional, and thus can not be well modeled by widely used dynamic models such as linear dynamical systems [8, 9] . Recently, many works [10, 11, 12, 13, 14, 15, 16] prefer to choose distributions of the gamma-Poisson family to build their hierarchical Bayesian models. In particular, these models enjoy strong explainability and can estimate uncertainty especially when the observations are noisy and incomplete. Among these works, Poisson-Gamma Dynamical Systems (PGDSs) [13] received a lot of attention because PGDS can learn how the latent dimensions excite each other to capture complicated dynamics in observed count series. For instance, a very inspiring research paper may motivate other researchers to publish papers on related topics [17] . The outbreak of COVID-19 in one state, may lead to the rapid rising of COVID-19 cases in the nearby states and vice versa [18] . In particular, PGDS can be efficiently learned with a tractable Gibbs sampling scheme via Poisson-Logarithmic data augmentation and marginalization technique [11] . Due to its strong flexibility, PGDS achieves better performance in predicting missing entities and future observations, compared with related models [9, 15] . Despite these advantages, PGDS still can not capture the time-varying transition dynamics underlying observed count sequences, which are commonly observed in real-world scenarios [19] . For instance, Submitted to 38th Conference on Neural Information Processing Systems (NeurIPS 2024). Do not distribute.during the initial stage of the COVID-19 pandemic, the worldwide counts of infectious patients were significantly affected by various local policies, government interventions, and emergent events [20, 21, 22] . The cross transition dynamics among the different monitoring areas were also evolving as the corresponding policies and interventions changed over time. Hence, PGDS unavoidably makes a certain amount of approximation error in capturing the aforementioned non-stationary count time series, using a time-invariant transition kernel.To mitigate this limitation, Non-Stationary Poisson-Gamma Dynamical Systems (NS-PGDSs), a novel kind of Poisson-gamma dynamical systems with non-stationary transition dynamics are developed.More specifically, NS-PGDS captures the evolving transition dynamics by the specifically-designed Dirichlet Markov chains. Via the Dirichlet-Multinomial-Beta data augmentation strategy, the Non-Stationary Poisson-Gamma Dynamical Systems can be inferred with a conjugate-yet-efficient Gibbs sampler. Our contributions are summarized as follows:• We propose a Non-Stationary Poisson-Gamma Dynamical System (NS-PGDS), a novel Poisson-gamma dynamical system with time-evolving transition matrices that can well capture non-stationary transition dynamics underlying observed count series.• Three Dirichlet Markov chains are dedicated to improving the flexibility and expressiveness of NS-PGDSs, for capturing the complex transition dynamics behind sequential count data.• Fully-conjugate-yet-efficient Gibbs samplers are developed via Dirichlet-Multinomial-Beta augmentation techniques to perform posterior simulation for the proposed Dirichlet Markov chains.• Extensive experiments are conducted on four real-world datasets, to evaluate the performance of the proposed NS-PGDS in predicting missing and future unseen observations. We also provide exploratory analysis to demonstrate the explainable latent structure inferred by the proposed NS-PGDS. 
preliminaries: 
 Let y (t) = y (t) 1 , • • • , y (t) V T ∈ N V be a vector of nonnegative count valued observations at time t.To capture the latent dynamics underlying count sequences, some previous works [23, 24] model the observations as y (t) = p z (t) , z (t) = f -1 x (t) , where p (•) is the observation likelihood function, and f (•) is an invertible link function that maps the parameters of observation component to continuous-valued latent variables x (t) ∈ R K . The latent factor x (t) evolves over time according to a linear dynamical system (LDS) given by x (t) ∼ N (Ax (t-1) , Λ -1 ), where A is the state transition matrix of size K ×K, andΛ = diag (λ 1 , • • • , λ K )is the inverse covariance matrix with λ -1 k determining the variance of k-th latent dimension. Han et al. [23] adopted the Extended Rank likelihood function to model count observations using LDS with time complexity O((K + V ) 3 ), which prevents it from practical applications for analyzing large-scale count data.Recently, Acharya et al. [15] and Schein et al. [13, 16] developed Poisson-gamma family models for sequential count observations. Gamma Process Dynamic Poisson Factor Analysis (GP-DPFA) [15] models count data as y(t) v ∼ Pois( K k=1 λ k ϕ vk θ (t)k ), where θ(t)k represents the strength of k-th latent factor at time t, and ϕ vk captures the involvement degree of k-th factor to v-th observed dimension.To ensure the model identifiability, we can impose a restriction as v ϕ vk = 1, and thus place a Dirichlet prior overϕ k = [ϕ 1k , • • • , ϕ V k ] T as ϕ k ∼ Dir (ϵ 0 , • • • , ϵ 0 ).To capture the underlying dynamics, the latent factor θ (t)k evolves over time according to a gamma Markov chain as θ(t) k ∼ Gam(θ (t-1) k , c t ), where c t is the rate parameter of the gamma distribution to control the variance of the gamma Markov chains. Although GP-DPFA can well fit one-dimensional count sequences, it fails to learn how the latent dimensions interact with each other.To address this concern, Schein et al. [13] developed Poisson-gamma dynamical systems to capture the underlying transition dynamics. In particular, θk evolves over time as θ(t) k ∼ Figure 1:The graphical representation of the NS-PGDS. The time interval is divided into equallyspaced sub-intervals. Each sub-interval contains M time steps. The transition dynamics is stationary within a sub-interval. In particular, the transition matrices evolve over sub-intervals via Dirichlet Markov processes while latent factors evolve over time steps via Eq. (1) .EQUATION, τ 0 ), where π kk2 represents how k 2 -th latent factor excites the k-th latent factor at next time step, andK k=1 π kk2 = 1.3 Non-Stationary Poisson-Gamma Dynamical Systems y (t) v ∼ Pois δ (t) K k=1 ϕ vk θ (t) k, in which, the latent factors are specified byEQUATIONwhere the multiplicative term δ (t) ∼ Gam (ϵ 0 , ϵ 0 ) and the transition matrices are time-varying asΠ (t) ≡ π (t) kk2 K k,k2=1. As shown in Figure 2 , to model the time-varying transition dynamics, we assume the whole time interval can be divided into I equally-spaced sub-intervals. The transition kernel behind complicated dynamic counts is assumed to be static within each sub-interval, while evolving over sub-intervals, to capture non-stationary behaviours. In another word, the proposed model allows the latent factors to evolve over time steps while the transition matrices change over sub-intervals but assumed to be stationary within each sub-interval, as shown in Figure 1 . In particular, we let each sub-interval contains M time steps, and the i-th interval contains time steps{t | t = (i -1) M + 1, • • • , iM }.We define i (t) as the function that maps time step t to its corresponding sub-interval.Dirichlet-Dirichlet Markov processes. To capture how the underlying transition kernel smoothly evolves over sub-intervals, we first propose the Dirichlet-Dirichlet (Dir-Dir) Markov chain asEQUATIONwhere π(i)k represents the k-th column of Π (i) , and the prior of the scaling parameter η is given by η ∼ Gam (e 0 , f 0 ).The initial states are defined as θ(1) k ∼ Gam (τ 0 ν k , τ 0 ). The prior for the transition kernel of the first sub-interval is given by π(1) k ∼ Dir (ν 1 ν k , • • • , ξν k , • • • , ν K ν k ), where ν k ∼ Gam( γ0 K , β) and ξ, β ∼ Gam (ϵ 0 , ϵ 0 ). Note that the expectation and variance of the transition kernel at i-th sub-interval can be calculated asE π (i) k | π (i-1) k = π (i-1) k , Var π (i) k1k | π (i-1) k = π (i-1) k1k 1 -π (i-1) k1k ηK + 1 ,respectively. The transition dynamics of i-th sub-interval inherits the information of the previous sub-interval, and also adapts to the data observed in the current sub-interval. The scaling parameter η controls the variance of the transition matrices.The prior specification defined in Eq.( 2) by rescaling the transition matrix at the previous sub-interval allows the transition dynamics to change smoothly, and thus might be insufficient to capture the rapid changes observed in complicated dynamics. To further improve the flexibility of the transition structure, two modified Dirichlet Markov chains are studied to capture the correlation structure between the dimensions of the transition matrices over time. Dirichlet-Gamma-Dirichlet Markov processes. We first introduce the Dirichlet-Gamma-Dirichlet (Dir-Gam-Dir)Markov chain to model the evolving transition matrices asπ (i) k ∼ Dir α (i) 1k , • • • , α (i) Kk , α (i) k1k ∼ Gam γ (i-1) k K k2=1 ψ (i-1) kk1k2 π (i-1)k2k , cEQUATIONwhere we use ψ(i-1)kk1k2 to capture the mutation between two consecutive sub-intervals, and its prior is given byEQUATIONEQUATIONk ∼ Gam (ϵ 0 , ϵ 0 ). Compared with the construction defined by Eq.( 2), the expectation of Dirichlet-Gamma-Dirichlet Markov chain isE π (i) k | π (i-1) k = Ψ (i-1) k π (i-1) k .This construction takes interactions among components of columns into account. Hence it will dramatically improve the flexibility of our model and thus better fit more complicated dynamics, compared with Dir-Dir Markov chains that only yield smoothing transition dynamics.Poisson-randomized-gamma-Dirichlet Markov processes. By leveraging the Poisson-randomized gamma distribution [25] , we introduce another type of time-varying transition kernels, which also model the interactions among components like Dir-Gam-Dir construction but may induce different properties such as sparsity. The Poisson-randomized-gamma-Dirichlet (PR-Gam-Dir) Markov chain can be formulated asEQUATIONwhere RG1 (•) denotes the randomized gamma distribution of the first type. Similarly, for ψEQUATIONk , and c(i)k , the priors are given byEQUATIONk ∼ Gam (ϵ 0 , ϵ 0 ) , respectively.The diagrams of three Dirichlet Markov constructions are shown in Figure 3 . 
theory assumptions and proofs: 
 Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes] Justification: The derivation of the Gibbs sampler is the main theoretical part of this paper which can be found in sec.4 and the appendix.Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced.• All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition.• Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.• Theorems and Lemmas that the proof relies upon should be properly referenced. 
markov chain monte carlo inference: 
 In this section, we present the Gibbs sampler for the proposed NS-PGDS. We only illustrate the key points of the derivation and the details can be found in the appendix.Lemma 1 If y ∼ NB (a, g (ζ)) and l ∼ CRT (y, a), where NB (•) refers to negative-binomial distribution, CRT (•) represents Chinese restaurant table distribution [26] , and g (z) = 1-exp (-z).Then the joint distribution of y and l can be equivalently distributed as y ∼ SumLog (l, g (ζ)) and l ∼ Pois (aζ) [11] , i.e.NB (y; a, g (ζ)) CRT (l; y, a) = SumLog (y; l, g (ζ)) Pois (l; aζ) ,where SumLog (l, g (ζ)) = l i=1 x i and x i ∼ Log (g (ζ)) are independently and identically logarithmic distributed random variables [27] .Lemma 2 Suppose n = (n 1 , • • • , n K ) and n | n ∼ DirMult (n, r 1 , • • • , r K ) ,where DirMult (•) refers to Dirichlet-multimonial distribution. We sample the augmented variableq | n ∼ Beta (n, r • ), where r • = K k=1 r k .According to [28] , conditioning on q, we haven k ∼ NB (r k , q).vk : Use the relationship between Poisson and multinomial distributions, we sampley (t) vk K k=1 | -∼ Mult   y (t) v , ϕ vk θ (t) k K k=1 ϕ vk θ (t) k K k=1   .Sampling ϕ k : Via Dirichlet-multinomial conjugacy, the posterior of ϕ k is(ϕ k | -) ∼ Dir ϵ 0 + T t=1 y (t) 1k , • • • , ϵ 0 + T t=1 y (t) V k .k : To sample from the posterior of θ l (t) k• | -∼ CRT y (t) •k + l (t+1) •k , τ 0 K k2=1 π i(t-1) kk2 θ (t-1) k2 , l (t) k1 , • • • , l (t) kK | -∼ Mult l (t) k• , π i(t-1) k1 θ (t-1) 1 K k2=1 π i(t-1) kk2 θ (t-1) k2 , • • • , π i(t-1) kK θ (t-1) K K k2=1 π i(t-1) kk2 θ (t-1) k2 . Let us define l (t) •k = K k1=1 l (t) k1k and ζ (t) = ln(1 + δ (t) τ0 + ζ (t+1)). After sampling the auxiliary variables, then for t = 1, • • • , T , by Poisson-gamma conjugacy, we obtainθ (t) k | -∼ Gam y (t) •k + l (t+1) •k + τ 0 K k2=1 π i(t-1) kk2 θ (t-1) k2 , τ 0 + δ (t) + ζ (t+1) τ 0 .Sampling Π (i) : We only illustrate Gibbs sampling algorithm for PR-Gam-Dir construction, sampling algorithms for other constructions can be found in the appendix. We define M as the length of each sub-interval, and I as the number of intervals.EQUATIONKk ) andEQUATION•Kk ) are multinomially distributed, where l(i) k1k = iM (i-1)M +1 l (t)k1k refers to the summation of l = 0, we sample the auxiliary variables as (qEQUATION•k ), then we have (l(i) k1k + g (i+1) •k1k ) ∼ NB(α (i) k1k , q (i) k ).Then we further sample (h(i) k1k | -) ∼ CRT(l (i) k1k + g (i+1) •k1k , α (i) k1k ). Via Lemma 1, we obtain h (i) k1k ∼ Pois(-α (i) k1k ln(1 -q (i) k )).For Dirichlet-Randomized-Gamma-Dirichlet Markov construction defined by Eq.( 4), we can equivalently represent it asEQUATIONWe define λ(i-1) k1k ≜ γ (i-1) k K k2=1 ψ (i-1) kk1k2 π (i-1) k2kfor notation conciseness. By Poisson-gamma conjugacy, we have (α(i) k1k | -) ∼ Gam(g (i) k1k + ϵ α + h (i) k1k , c (i) k -ln(1 -q (i) k )). If ϵ α > 0, we can sample the posterior of g (i) k1k via (g (i) k1k | -) ∼ Bessel(ϵ α -1, 2 α (i) k1k c (i) k λ (i-1)k1k ), where Bessel (•)denotes Bessel distribution. If ϵ α = 0, we sample g (i) k1k via g (i) k1k | -∼        Pois c (i) k λ (i-1) k 1 k c (i) k -ln 1-q (i) k if h (i) k1k = 0 SCH h (i) k1k , c (i) k λ (i-1) k 1 k c (i) k -ln 1-q (i) k otherwise,where SCH (•) denotes the shifted confluent hypergeometric distribution [16] . Defining g(i) k1k = g (i) k1•k = K k2=1 g (i)k1k2k , we first augmentg (i) k11k , • • • , g (i) k1Kk ∼ Mult g (i)k1k , ψ(i-1) kk1k2 π (i-1) k2k K k2=1, then we obtain gEQUATION). By Dirichlet-multinomial conjugacy, we haveψ (i-1) k1k2 , • • • , ψ (i-1) kKk2 | -∼ Dir ϵ 0 + g (i) 1k2k , • • • , ϵ 0 + g (i)Kk2k , andπ (i-1) k | -∼ Dir α (i-1) 1k + l (i-1) 1k + g (i) •1k , • • • , α (i-1) Kk + l (i-1) Kk + g (i)•Kk .Specifically, we have α  and α (1)(1) k1k = ν k1 ν k , if k 1 ̸ = k,k1k = ξν k , if k 1 = k.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes] Justification: We carefully described the proposed model in sec.2 and sec.3 and the experiment details can be found in sec.6.1.Guidelines:• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.• If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.• Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.• While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 
related work: 
 Modeling count time sequences has been receiving increasing attentions in statistical and machine learning communities. Han et al. [23] adopted linear dynamical systems to capture the underlying dynamics of the data and leveraged Extended Rank likelihood function to model count observations. Some Poisson-gamma models assume that the count vector at each time step is modeled by Poisson factor analysis (PFA) [11] and leverage special stochastic processes to model the temporal dependencies of latent factors. For example, gamma process dynamic Poisson factor analysis (GP-DPFA) [15] adopts gamma Markov chains which assumes the latent factor of the next time step is drawn from a gamma distribution with the shape parameter be the latent factor of the current time step. Schein et al. [13] proposed Poisson-gamma dynamical systems (PGDSs), which take the interactions among latent dimensions into account and use a transition matrix to capture the interactions. Deep dynamicPoisson factor analysis (DDPFA) [29] adopts recurrent neural networks (RNNs) to capture the complex long-term dependencies of latent factors. Yang and Koeppl [30] applied Poisson-gamma count model to analyze relational data arising from longitudinal networks, which can capture the evolution of individual node-group memberships over time. Many modifications of PGDS have been proposed in recent years. Guo et al. [31] proposed deep Poisson-gamma dynamical systems which aim to capture the long-range temporal dependencies. Schein et al. [16] employed Poisson-randomized gamma distribution to build a new transition process of latent factors. Chen et al. [32] proposed Switching Poisson-gamma dynamical systems (SPGDS), allowing PGDS to select from several transition matrices, and thus can better adapt to nonlinear dynamics. In contrast to SPGDS, the number of transition matrices of the proposed NS-PGDS is not limited and thus can be adopted to analyze various complicated non-stationary count sequences. Filstroff et al. [33] extensively analyzed many gamma Markov chains for non-negative matrix factorization and introduced new gamma Markov chains with well-defined stationary distribution (BGAR).Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [No] Justification: The authors will release the data and code as soon as possible if this paper could be accepted.Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details.• While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).• The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: //nips.cc/public/guides/CodeSubmissionPolicy) for more details.• The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.• The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why.• At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).• Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
experiments: 
 We conducted experiments for both predictive and exploratory analysis to demonstrate the ability of the proposed model in capturing non-stationary count time sequences. The baseline models includedin the experiments are: 1) Gamma process dynamic Poisson factor analysis (GP-DPFA) [15] .GP-DPFA models the evolution of latent components as θ ).(t) k ∼ Gam(θ (t-1) k , c t ),3) Gamma Markov chains on the rate parameter with hierarchical auxiliary variable (GMC-HIER) [33] . GMC-HIER models the evolution of latent components with an auxiliary variables as z(t) k ∼ Gam(α z , β z θ (t-1) k) and θ(t) k ∼ Gam(a θ , β θ z (t)k ). 4) Autogressive beta-gamma procecss (BGAR) [34, 33] . BGAR is also a gamma Markov model. In contrast to the above models, there is a well-defined stationary distribution for BGAR. 5) Poisson-gamma dynamical system (PGDS) [13] takes interactions among latent dimensions into account, and models the evolution of θ(t) k as θ (t) k ∼ Gam(τ 0 K k2=1 π kk2 θ (t-1) k2, τ 0 ).The real-world datasets used in the experiments are: 1) Integrated Crisis Early Warning System (ICEWS): ICEWS is an international relations event dataset, comprising interaction events between countries extracted from news corpora. For ICEWS dataset, we have T = 365 time steps and V =To compare the predictive performance of the proposed model with the baselines, we considered two standard tasks: data smoothing and forecasting. For data smoothing task, our objective is to predicty (t)given the remaining data observation Y \y (t) . To this end, we randomly masked 10 percents of the observed data over non-adjacent time steps, and predicted the masked values. For forecasting task, we held out data of the last S time steps, and predicted y (T +1) , • • • , y (T +S) given y (1) , • • • , y (T ) . In this experiment we set S = 2. We ran the baseline models including GP-DPFA, PGDS, GMC-RATE, GMC-HIER, BGAR, using their default settings as provided in [15, 13, 33] . For the NS-PGDS, we set K = 100 for ICEWS, K = 10 for other datasets, and set τ 0 = 1, γ 0 = 50, ϵ 0 = 0.1. We performed 4000 Gibbs sampling iterations. In the experiments, we found that the Gibbs sampler started to converge after 1000 iterations, and thus we set the burn-in time be 2000 iterations. We retained every hundredth sample, and averaged the predictions over the samples. Mean relative error (MRE) and mean absolute error (MAE) are adopted to evaluate the model's predictive capability, which are defined as MRE = 1T V t v |y (t) v -ŷ (t) v | 1+y (t) v and MAE = 1 T V t v | y (t) v - ŷ(t) v | respectively,where y(t)v indicates the true count and ŷ(t) v is the prediction.As the experiment results shown in Table 1 , the NS-PGDS exhibits improved performance in both data smoothing and forecasting tasks. We attribute this enhanced capability to the time-varying transition kernels, which effectively adapt to the non-stationary environment, and thus achieve improved predictive performance. For some datasets (e.g. ICEWS) and tasks, the effectiveness of the Dir-Gam-Dir and Pr-Gam-Dir constructions does not be exhibited in the numerical results. However, these two constructions indeed induce more informative patterns compared with Dir-Dir construction, as shown in the exploratory analysis. "network-neural-networks", "kernelimage-space", "network-neuralnetworks", "neural-networks-state".The representative words for latent factors in the second block are "image-sparse-matrix", "kernelsupervised-random", "matrix-sample-random", "inference-prior-latent", "state-policy-gamma". The first block primarily captured the correlations among the research topics about neural networks.The second block reflects that, from the 1990s, statistical learning and Bayesian methods began to dominate, and these topics are highly correlated. Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: The experiment details can be found in sec.6.1.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.• The full details can be provided either with the code, in appendix, or as supplemental material. 
conclusion: 
 The Poisson-gamma dynamical systems with time-varying transition matrices, have been proposed to capture complicated dynamics observed in non-stationary count sequences. In particular, Dirichlet Markov chains are constructed to allow the underlying transition matrices to evolve over time.Although the Dirichlet Markov processes lack conjugacy, we have developed tractable-but-efficient Gibbs sampling algorithms to perform posterior simulation. The experiment results demonstrate the improved performance of the proposed NS-PGDS in data smoothing and forecasting tasks, compared with the PGDS with a stationary transition kernel. Moreover, the experimental results on several real-world data sets show the explainable structures inferred by the proposed NS-PGDS. For the future work, we plan to design a method that can find the point of change and thus the length of each sub-interval can be determined automatically instead of a constant. We also consider to generalize Dirichlet belief networks by incorporating the proposed Dirichlet Markov chain constructions, which allow the hierarchical topics to mutate across layers, and thus can generate more rich text information.And we also consider to capture non-stationary interaction dynamics among individuals over online social networks in the future research.Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes]Justification: Table 1 reports the predictive performance of the proposed model and the corresponding standard deviation. The results are computed by running the Gibbs sampling several times from different initialization.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).• The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)• The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean.• It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified.• For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).• If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: The experiments are conducted on a server with an Intel(R) Xeon(R) CPU E5-2699Cv4 @ 2.20GHz and 64G RAM.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.• The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.• The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.• The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes] Justification: For positive societal impacts, we have discussed in conclusion section for the potential application for textual analysis and social networks. And the authors think this work does not have potential negative societal impacts.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.• The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.• If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).• At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification: This paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.• Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper.• Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.• Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.• We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.• For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: Notation. When expressing the full conditionals for Gibbs sampling, we use the shorthand "-" to denote all other variables. We use "•" as an index summation shorthand, e.g., x •j = i x ij .In this section, we present a fully-conjugate and efficient Gibbs sampler for the proposed NS-PGDS.The sampling algorithms depend on several key technical results, which we will repeatedly exploit, thus we list them below. Lemma 2. Suppose n = (n 1 , • • • , n K ) and n | n ∼ DirMult (n, r 1 , • • • , r K ) , where DirMult (•)refers to Dirichlet-multimonial distribution. We sample the augmented variable q | n ∼ Beta (n, r • ),where r • = K k=1 r k . According to [28] , conditioning on q, we have n k ∼ NB (r k , q). (y 1 , • • • , y S ) ∼ Mult(y • , ( λ1 S s=1 λs , • • • , λ S S s=1 λs )), where Mult (•) represents multinomial distribution [35] .vk : Use the relationship between Poisson and multinomial distributions as described by Lemma 3, given observed counts and latent parameters, we sampleEQUATIONThen the distribution of y(t)vk is y(t) vk ∼ Pois(δ (t) ϕ vk θ (t) k ).Sampling ϕ k : Via Dirichlet-multinomial conjugacy, the posterior of ϕ k isEQUATIONMarginalizing over θ(t) k : Note that y (t) v = y (t) v• = K k=1 y (t)vk and y(t) vk ∼ Pois(δ (t) ϕ vk θ (t) k ). Then we define y (t) •k = V v=1 y (t) vk . Because V v=1 ϕ vk = 1, we obtain y (t) •k ∼ Pois(δ (t) θ (t) k ).We start by marginalizing over θ(T )k , using the definition of negative-binomial distribution, we obtainy (T ) •k ∼ NB τ 0 K k2=1 π i(T -1) kk2 θ (T -1) k2 , g ζ (T ) ,where ζ (T ) = ln(1 + δ (T ) τ0 ). Next, we further marginalize over θ (T -1) k. To this end, we first sample auxiliary variablesEQUATION.By Lemma 1, the joint distribution of y (T )•k and l (T ) k can be expressed asEQUATION.Via Lemma 3, we re-express the auxiliary variables asl (T ) k = l (T ) k• = K k2=1 l (T )kk2 , and obtain l(T ) kk2 ∼ Pois ζ (T ) τ 0 π i(T -1) kk2 θ (T -1) k2.Then we define l (T )•k = K k1=1 l (T )k1k . Leveraging Lemma 3 andK k1=1 π i(T -1) k1k = 1, we obtain l (T ) •k ∼ Pois ζ (T ) τ 0 θ (T -1) k and l (T ) 1k , • • • , l (T ) Kk ∼ Mult l (T ) •k , π i(T -1) 1k , • • • , π i(T -1) Kk .Next, note that y(T -1) •k ∼ Pois(δ (T -1) θ (T -1) k ), if we introduce m (T -1) k = y (T -1) •k + l (T )•k , then we havem (T -1) k ∼ Pois θ (T -1) k δ (T -1) + ζ (T ) τ 0 .Because the prior of θ(T -1) kis gamma distributed, by the definition of negative-binomial distribution, we can again marginalize over θ(T -1) k to obtain m (T -1) k ∼ NB τ 0 K k2=1 π i(T -2) kk2 θ (T -2) k2 , g ζ (T -1),whereζ (T -1) = ln(1 + δ (T -1) τ0 + ζ (T )). Then we introduce auxiliary variablesl (T -1) k ∼ CRT m (T -1) k , τ 0 K k2=1 π i(T -2) kk2 θ (T -2) k2.And similar to the case for t = T , we can obtainl (T -1) •k ∼ Pois ζ (T -1) τ 0 θ (T -2) k and m (T -2) k ∼ NB τ 0 K k2=1 π i(T -3) kk2 θ (T -3) k2 , g ζ (T -2) .Thus we have marginalized over θ EQUATIONl (t) k1 , • • • , l (t) kK | -∼ Mult l (t) k• , π i(t-1) k1 θ (t-1) 1 K k2=1 π i(t-1) kk2 θ (t-1) k2 , • • • , π i(t-1) kK θ (t-1) K K k2=1 π i(t-1) kk2 θ (t-1) k2 . (9)And via Lemma 3, we obtain -1 ) Kk (10) We compute ζ (t) recursively vial (t) 1k , • • • , l (t) Kk ∼ Mult l (t) •k , π i(t-1) 1k , • • • , π i(tEQUATION)After sampling the auxiliary variables, then for t = 1, • • • , T , by Poisson-gamma conjugacy, we obtainEQUATION•k + l(2)EQUATIONEQUATION)Sampling Π (i) : We define M as the length of each sub-interval, and I as the number of intervals.For i = I, by Eq.( 10), (l(I) 1k , • • • , l (I)Kk ) is multinomial distributed. Thus by multinomial-Dirichlet conjugacy, we obtainEQUATIONwhere l Inference for Dirichlet-Dirichlet Markov chains. For Dirichlet-Dirichlet Markov chains, α(i) k1k = ηKπ (i-1)k1k . By Eq.( 10), (l(i) 1k , • • • , l (i) Kk ) is multinomial distributed. If we marginalize (π (i) 1k , • • • , π (i) Kk ), (l (i) 1k , • • • , l (i)Kk ) will be Dirichlet-multinomial distributed. Thus by Lemma 2, for i = I, we first sample the auxiliary variables asq (I) k | -∼ Beta l (I)•k , ηK and h(I) k1k | -∼ CRT l (I)k1k , ηKπEQUATION)Similarly, by Eq.( 18), (h(i) 1k , • • • , h (i)Kk ) is also Dirichlet-multinomial distributed. Thus for i =I -1, • • • , 2,we sample the auxiliary variables asq (i) k | -∼ Beta l (i) •k + h (i+1) •k, ηK and h(i) k1k | -∼ CRT l (i) k1k + h (i+1)k1k , ηKπEQUATIONwhere l(i) k1k = iM (i-1)M +1 l (t)k1k refers to the summation of l(t)k1k over i-th interval. Via Lemma 2, conditioning on q(i)k , we haveEQUATION.Then via Lemma 1, we obtainEQUATION)Note that by Eq.( 17), h(i)k1k is Poisson distributed and by Lemma 3, we obtainEQUATION)In addition, note thatEQUATION)Specifically, for i = 1, we haveEQUATION)For sampling η, note that (h(i) k1k | -) ∼ Pois(-ηKπ (i-1) k1k ln 1 -q (i) k ), i = I, • • • , 2.Given the prior η ∼ Gam (e 0 , f 0 ), via Poisson-gamma conjugacy, we obtainEQUATION)Inference for Dirichlet-Gamma-Dirichlet Markov chains. For Dirichlet-Gamma-Dirichlet Markov chainsEQUATION.By Eq.( 10), (l(i) 1k , • • • , l (i) Kk ) is multinomial distributed. If we marginalize (π (i) 1k , • • • , π (i) Kk ), (l (i) 1k , • • • , l (i)Kk ) will be Dirichlet-multinomial distributed. Thus by Lemma 2, for i = I, we first sample the auxiliary variables asEQUATIONSimilarly, by Eq.( 27), (gEQUATION•Kk ) is also Dirichlet-multinomial distributed. Thus for i =I -1, • • • , 2,we sample the auxiliary variables asEQUATION)Via Lemma 2, conditioning on q(i)k , we haveEQUATION.Then via Lemma 1, we obtainh (i) k1k ∼ Pois -α (i) k1k ln 1 -q (i) k.Thus via Poisson-gamma conjugacy, we obtainEQUATION)Marginalizing over α(i)k1k , and via the definition of negative-binomial distribution, we haveEQUATIONThen using Lemma 1, we sampleEQUATIONand obtaing (i) k1k ∼ Pois γ (i-1) k K k2=1 ψ (i-1) kk1k2 π (i-1) k2k ln 1 -ln 1 -q (i) k c (i) k . If we define g (i) k1k = g (i) k1•k = K k2=1 g (i)k1k2k , and augmentEQUATION.By Lemma 3, we haveg (i) k1k2k ∼ Pois γ (i-1) ψ (i-1) kk1k2 π (i-1) k2k ln 1 -ln 1 -q (i) k c (i) k .Using Lemma 3 andEQUATION•Kk ∼ Mult gEQUATION(i) 1k2k , • • • , g (i) Kk2k ∼ Mult g (i) •k2k , ψ (i-1) kk1k2 K k1=1 .Thus by Dirichlet-multinomial conjugacy, for i = I, • • • , 2, we can obtainEQUATIONEQUATIONFor sampling γ(i-1) k, note that by Eq.( 26) andK k1 ψ (i-1) kk1k2 = 1, we have g (i) •k = K k1=1 g (i)k1k and gEQUATIONThus via Poisson-gamma conjugacy, we obtainEQUATIONBy gamma-gamma conjugacy, we haveEQUATIONInference for Dirichlet-Randomized-Gamma-Dirichlet Markov chains.For Dirichlet-Randomized-Gamma-Dirichlet Markov chains,EQUATION, which can be equivalently represented asEQUATION, and gEQUATION.By Eq.(Kk ) is multinomial distributed. If we marginalize (πEQUATIONKk ) will be Dirichlet-multinomial distributed. Thus by Lemma 2, for i = I, we first sample the auxiliary variables asEQUATION•k and hEQUATIONSimilarly, by Eq.(39), (gEQUATION•Kk ) is also Dirichlet-multinomial distributed. Thus for i =I -1, • • • , 2,we sample the auxiliary variables asEQUATION)Via Lemma 2, conditioning on q (i) k , we haveEQUATION.Then via Lemma 1, we obtainh (i) k1k ∼ Pois -α (i) k1k ln 1 -q (i) k .Via Poisson-gamma conjugacy, we first sampleEQUATION)If ϵ α > 0, we can sample the posterior of gEQUATIONwhere Bessel (•) denotes Bessel distribution. If ϵ α = 0, we sample gEQUATIONJustification: As we discussed in the conclusion part, the length of each sub-interval is a constant and is treated as a hyper-parameter of the model. In the future work, we plan to design a method that can find the point of change and thus the length of each sub-interval can be determined automatically.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.• The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.• The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.• The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.• The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.• If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.• While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 


Decision: The reviewers agree that this is a fairly solid applied paper, solving a natural problem: how can an existing model PGDS be extended to nonstationary settings. However, the proposed approach is a bit incremental (piecewise-constant nonstationarity), with the choice of how these piecewise intervals are chosen a bit ad hoc. There were also concerns about the experimental evaluation. Overall, this paper falls a bit below the threshold set by some of the more competitive submissions to NeurIPS.

Reviews: review: paper_summary: paper_summary: This paper introduces non-stationary Poisson-Gamma dynamical systems, an extension of Poisson Gamma dynamical systems with a dynamic transition matrix. Decomposing the time steps into equally spaced subintervals, the transition matrices evolve between sub-intervals, remaining static within sub-intervals. The authors introduce three options for transitions to occur. The authors derive a Gibbs sampling scheme for exact posterior inference using data augmentation techniques and showcase the effectiveness of their method through a series of predictive and qualitative results.
strength_weakness: strength_weakness: I'm not convinced that the magnitude of the author's contribution, nor the significance of the paper is strong enough to warrant acceptance, and the methods produce only marginally better results than that of Poisson Gamma dynamical systems. The qualitative results are not groundbreaking.
questions: questions: How does inference scale with the number of sub-intervals in each of the three methods?
How does the user choose the number of sub-intervals? Can they be fit adaptively?
How much extra training time do the proposed methods add, as compared to Poisson Gamma dynamical systems?
limitations: limitations: The authors address the limitations of their work, stating intention to address these limitations (e.g. constant sub-interval lengths) in future work.
paper_summary: paper_summary: The work extends  Poisson-Gamma Dynamical systems (PGDS) to model non-stationary dynamics by replacing the constant transition matrix $\Pi$ with a time dependent one $\Pi^{(t)}$ and the original Dirichlet prior on the columns with three different Dirichlet Markov chain constructions. The manuscript describes am efficient Gibbs-sampler for inference.
strength_weakness: strength_weakness: My major problem is the experiment evaluation. In Table 1 in the NIPS dataset we can see results like $14.014 \pm 4.387$ bolded, over values like $14.706 \pm 4.414$, or $17.105 \pm 6.449$. 
In ICEWS values like $0.214 \pm 0.008$ over $0.215 \pm 0.007$ , in USEI $4.596 \pm 0.562$ over $4.703  \pm 0.538$, in COVID $6.969 \pm 1.107$ over $7.566 \pm 1.095$. These are mainly smoothing results. In the light of this I am not confident in the statement “As the experiment results shown in Table 1, the NS-PGDS exhibits improved performance in both data smoothing and forecasting tasks.”.  We do not know how the confidence interval was computed, or how many repeats were made. The lack of statistical rigor in the evaluation stands in striking contrast with the sophisticated Bayesian model presented. 

Besides this, other possible problem with the evaluation is that the manuscript states that default paramerters were used for the benchmark methods “GP-DPFA, PGDS, GMC-RATE, GMC-HIER, BGAR”  while the present method used specific K based on the dataset. It is very hard to tell if this is a fair comparison or not.
questions: questions: 1)  The time series are divided  to equally-spaced sub-intervals, and a $\Pi$ transition matrix is infered for all interval. Why not use a diferent matrix for all time points with a very slowly varying Markow chain. Is this decision was made due to computational reasons or due to modeling assumptions?
2)  Please provide the parameters for the benchmark methods in the appendix to facilitate comparison, it is very time consuming for the reader to search all default parameters from the references.
3) Please bold all indistinguishable results in Table 2 using a statistical test, and reassess your conclusions.
4) What is the reason of duplicated  factors in the exploratory analysis like “neural-network-networks” and “network-neural-networks”. What the order of the words means? Strength of association?
limitations: limitations: No specific limitation section was provided. The part on future work in the Conclusion can be interpreted as pointing out some limitations of the current model, but a specific limitation statement would be preferable.
paper_summary: paper_summary: Existing PGDS models struggle with capturing the time-varying transition dynamics seen in real-world data. To address this, the submission proposed a non-stationary PGDS, allowing the transition matrices to evolve over time, modeled by Dirichlet Markov chains. Using Dirichlet-Multinomial-Beta data augmentation techniques, a fully-conjugate and efficient Gibbs sampler is developed for posterior simulation. Experiments demonstrate that the proposed non-stationary PGDS achieves improved predictive performance compared to related models.
strength_weakness: strength_weakness: I did not find any obvious weaknesses.
questions: questions: None.
limitations: limitations: The authors discussed some future work directions in the conclusion.
paper_summary: paper_summary: This work extends Poisson-Gamma Dynamical Systems (PGDSs) by considering non-stationary transition dynamics to effectively capture the evolving dynamics of observed count sequences.

The authors propose a model where the underlying transition matrices evolve over time, based on three (gradually more complex and flexible) Dirichlet Markov chains.

For inference of the model, the authors make use of the Dirichlet-Multinomial-Beta data augmentation to derive a fully-conjugate Gibbs sampler.

Experiments showcase improved data-smoothing and forecasting performance of the proposed method across several real-world datasets.
strength_weakness: strength_weakness: - The main limitation of this work is the assumption that the transition kernel is static within each sub-interval: i.e., the authors consider that the kernel can only change at discrete instants, while is constant within each sub-interval.
questions: questions: - Can the authors justify and explain their choice for only allowing discrete-time transition kernel changes?
    - How can one determine the length of sub-intervals in practice? How did the authors determine these sub-intervals in their experiments?
    - Would it be possible to accommodate sub-intervals of varying length, $M$, and what would be the implications?
    - Would it be possible to consider a continuously changing transition kernel? What would be the implications for the model and/or the estimation procedure?

- The different transition kernel evolution models proposed do not only differ in their flexibility to capture different phenomena, but also on their complexity:
    - Can the authors elaborate on the number of learnable parameters of each model?
    - What is the computational and statistical complexity associated with each?
    - Results do not seem to provide data-smoothing and forecasting performance improvements: is the added flexibility worth the complexity?

- Sections 5 seems to have quite an overlap with some of the preliminaries introduced in Section 2:
   - Would it be possible to merge both, or is there a reason why these two should be self-contained in different sections?
limitations: limitations: The authors address the main limitations of their work.

Comments: title: Thank you for the clarification
Thank you!
comment: Dear Reviewer:

Sorry for misunderstanding your comments, and thanks for your response. Compared with the forecasting task, data smoothing is much simpler and PGDS is competent to this task very well. More importantly, for data smoothing task, we randomly masked 10 percents of the observed data over non-adjacent time steps, and predicted the masked values. The random selection of masked data also poses significant variance of the numerical results. Therefore, it is **very hard** for NS-PGDS to outperform PGDS in data smoothing task by a large margin.

Besides,  the main motivation of this work is to allow the transition matrix of PGDS to be time-varying. The results of forecasting task and exploratory analysis have partially validated the effectiveness of the proposed method. Though for data smoothing, the numerical results are not that significant in statistical, however, at least, the results of NS-PGDS are not worse than that of PGDS, and this does not conflict with our motivation.

As the reviewer has pointed out, we will reassess our conclusions and give more detailed explanations about results of data smoothing task in the final version.
Thank you for your response. I believe these details add to the quality of the contribution, although it should be noted that the smoothing results are not significant! The superior forecasting results, however, are convincing. The little extra training time and superior forecasting results increase my confidence in this work, and I will change my recommendation accordingly.
Thank you for the detailed response! My point was not that the method outperforms in only smoothing result, quite the opposite: that the listed questionable claims (that you yourself shown by the p-value analysis to be insignificant) i mentioned are smoothing results.
I thank the reviewers for their response to my questions, specifically, on the complexity details of each algorithm.