Bayesian quantile regression for streaming data

1.   Introduction

The term "big data" may have different meanings in different fields, so it has thus become a dominant topic in almost all academic and application fields. In a more general sense, big data encapsulates aspects such as quantity, diversity, speed, variability, and accuracy ^[1]. However, applying statistical models and methods to big data often entails a substantial computational burden, straining computer memory due to its large size and impacting computational efficiency, in that even seemingly simple tasks can demand considerable computation time. In order to solve the dilemma of big data analysis, some researchers have proposed three methods ^[2]: Subsampling-based approaches ^[3,4], divide and conquer approaches ^[5,6,7,8], and online updating approaches ^{[9,10,11,12,13,14]}. Differing from the other two methods, online updating approaches focus on streaming data and address statistical problems in an updating framework without storage requirements for historical data. In the era of big data, streaming datasets have gradually appeared in various fields, including investment analysis, medical imaging, and computervision, making streaming data analysis increasingly useful and vital across various industries. For instance, financial institutions track the changes of the stock market in realtime, calculate the value at risk, and automatically rebalance the portfolio according to the stock price movements. Due to the increasing demand for stream processing, online updating became more particularly appealing with its ability to quickly process huge volumes of data so that organisa tions or businesses can react to changing conditions in realtime.

In the framework for regression-type analysis, Shi et al. ^[11] developed online updating algorithms for linear models and estimating equations. However, theoretical gurantees of these methods have been established based on some strong regularity conditions. Mohamad and Bouchachia ^[13] proposed a renewable estimation for the generalized linear model, which overcame the aforementioned limitations. Lin et al. ^[15] introduced a general framework to execute renewable weighted sum estimation in various online updates. Wang et al. ^[16] expanded the scope of online update methods by adapting to new predictive variables midway through the data flow. Wu and Chen ^[17] developed an online update method for survival analysis under the Cox proportional hazard model. Furthermore, Xue et al. ^[18] proposed a online updating approach to evaluate the proportional hazard hypothesis. Balakrishnan and Madigan ^[19] proposed a one-pass streaming algorithm for Bayesian regression. Based on sketch drawing technology and random projection, Geppert et al. ^[20] improved the streaming algorithm for Bayesian regression proposed in ^[19]. All of the above methods and algorithms are applicable in practice. However, they mainly focus on mean regression, or are built upon likelihood framework. Moreover, it is noteworthy that mean regression and likelihood methods still have limitations in robustness, and are very sensitive to outliers or heavy-tailed distributions.

As a natural upgrade solution, quantile regression (QR) ^[21] has been used as a feasible alternative to mean regression and has become a popular data analysis strategy. Quantile regression does not simply focus on the average value of the distribution but considers the entire conditional distribution, thus being more comprehensive. In addition, it is also robust to heavy-tailed distributions and outliers. Because of these advantages, quantile regression has a wide range of applications, namely in health research ^[22], longitudinal data analysis ^[23,24], economics ^[25], and machine learning ^[26,27]. However, because the datasets arrive continuously in an unbounded stream and the computer memory is usually limited, the traditional quantile regression method cannot retain them when processing stream data. Therefore, the traditional quantile regression strategy, which assumes that the computer memory can adapt to all of the datasets, is no longer applicable. In addition, the estimator obtained by quantile regression has no closed form, and the quantile regression loss function is also nondifferentiable. This also leads to the fact that streaming data analysis methods cannot be directly applied to quantile regression; see Chen et al. ^[28].

To circumvent the nondifferentiability of the QR loss function, Chen et al. ^[28] introduced Horowitz's smoothed QR, involving smoothing the indicator part of the check function through a kernel smoothing survival function. Additionally, Chen and Zhou ^[29] proposed online QR strategies, employing a smooth function to approximate the nondifferentiable check loss function. The method in ^[28] introduced a linear-type estimator for QR, and its calculation process needs to select bandwidth parameters through data-driven criteria. Furthermore, the approach in ^[29] relied on the Taylor expansion of the estimating equations, requiring the estimation of the unknown conditional density function of the response variable. More recently, Wang et al. ^[30] developed a renewable QR estimation strategy for streaming datasets, which renews the estimator with current data and summarized statistics of historical data rather than historical raw data. In summary, the above methods have effectively solved the dilemma of quantile regression in the case of streaming data, but these methods focus on the solution of frequentist rather than Bayesian. In the field of Bayesian streaming data, current research work is mainly focused on Bayesian flow networks (BFN)^[48,49]. Rezende and Mohamed ^[48] introduced the concept of normalizing flows, which is a key component of Bayesian flow networks. Müller and Quintana^[49] provided a comprehensive introduction to the application of Bayesian nonparametric models in data analysis and inference, including streaming data. BFN is a generative model that integrates the advantages of Bayesian inference and deep learning. It has received extensive attention in recent years in areas such as data generation, image generation, and text generation. Bayesian flow networks particularly excel in handling discrete data generation problems, as they can generate continuous and differentiable samples, providing a new solution for discrete data generation. BFN combines the streaming models of deep learning with Bayesian inference, enabling the model to adapt to both continuous and discrete data, and enhancing the interpretability and flexibility of generative models. In tasks such as image generation and text generation, BFN has achieved good experimental results, demonstrating its tremendous potential in practical applications. Therefore, it is interesting to explore a more convenient Bayesian quantile regression strategy for streaming datasets. To the best of our knowledge, there is little related work on Bayesian quantile regression with streaming data.

In this paper, in contrast to the above methods, we develop a new idea for the analysis of streaming data using the Bayesian quantile model, which renews the posterior distribution with current data and summary statistics of historical data, instead of historical raw data. Specifically, we only need to retain the posterior distribution parameters of previous data and a likehood function for current data stream, and the historical data can be completely discarded. The Bayesian streaming quantile regression (BSQR) model is formulated for streaming data through a conjugate normal-inverse-gamma (NIG) prior distribution. Leveraging the posterior distribution of historical data as the prior for current data stream, a posterior distribution is established without any loss of information. Theoretically, we prove that the Bayesian streaming quantile regression estimator is equivalent to the global oracle estimator calculated based on the entire data. Numerical experiments on both synthetic and real data are also included to validate the theoretical results and illustrate the good performance of our new method.

The rest of the paper is organized as follows. Section 2 presents the introduction of BQR method. Section 3 gives the detailed expressions of NIG prior and posterior distributions for BQR using an informative $g$-prior. The BSQR method is developed in Section 4. Sections 5 and 6 validate the obtained theoretical results and illustrate the good performance of the proposed BSQR approach on synthetic and real data. Section 7 contains the discussion and conclusions. In addition, the proof of Theorem 1 is provided in the appendix.

2.   Bayesian quantile regression model

2.1.   Quantile regression based on AL distribution

We consider the following linear quantile regression model

where $y_i$ is a continuous response variable and $\mathit{\boldsymbol{x_i}}$ is an $k \times 1$ vector of predictors for the $i$-th observation, $i = 1, \cdots, n$. $\beta$ is an $k \times 1$ vector of unknown parameter of interest, $e_i$ is the error term, and its distribution is assumed to have zero $\tau$-th quantile.

Now, consider a linear conditional quantile function as follows:

where $Q_{\tau}(\cdot) = inf\{y:F_{y_i}(y|x_i)\ge \tau\}$ denotes the inverse of the cumulative distribution function of the response conditional on $x_i$. In $\tau$-th quantile regression, the estimation for $\beta$ can be obtained by minimizing

where $\rho_{\tau}(u) = u(\tau-I(u\leq 0))$ denotes the check function.

Koenker and Machado ^[50] showed that $\rho_{\tau}(u)$ is exactly matched to the asymmetric Laplace distribution (ALD), whose density function is given by

where $\mu$ is the location parameter and $\sigma$ is the scale parameter.

Assume that errors $\{e_i\}_{i = 1}^n$ are i.i.d $ALD(0, \sigma, \tau)$, the likehood of $y_i$ is

According to Koenker and Machado ^[50], minimizing (2.3) is equivalent to maximizing the above likelihood function (2.5). Yu and Moyeed ^[51] proposed a Bayesian method for QR by assuming an ALD on the response as a working likelihood. Over the past few years, this working likelihood function has been used in Geraci and Geraci and Bottai ^[52,53]; among others. So, it is a natural choice to use the ALD error distribution in Bayesian quantile models; "we believe that the ALD approach is a valuable and relatively simple alternative to these other methods that does not need complex choices of prior distributions and prior parameters"; see Section 3 of Benoit et al.^[54].

However, direct use of the likelihood function is rather inconvenient for Bayesian quantile regression. Specifically, the conditional distribution for the regression coefficients is not analytically tractable due to the complexity of the above likelihood function. By Chu et al. ^[31] and Lum et al. ^[32], the $e_i$ can be denoted as a location-scale mixture of normals as follows:

where, $z_i \sim Exp(\sigma^{-1})$, $u_i \sim N(0, 1)$, and $z_i$ and $u_i$ are independent. Let $v_i = \frac{1}{\tau(1-\tau)}z_i$, then

Then, we have

where, $Y = \{y_1, y_2, \cdots, y_n\}^T$ is an $n \times 1$ response vector of $y_i$, $X$ is an $n \times k$ predictor matrix with ${\it i}$-th row $\mathit{\boldsymbol{x}}_i^T$, $V$ denotes an $n \times 1$ vector of $v_i$, and $\Sigma$ is the diagonal matrix of $v_i$.

Subsequently, let $Y_{\tau}^{*} = \frac{1}{\sqrt{2}}\left(Y-\left(1-2\tau \right)V \right)$ and ${{X}^{*}} = \frac{1}{\sqrt{2}}X$, we have

and

2.2.   Likelihood function and prior distribution

According to Chu et al. ^[31], a NIG distribution is given by the following definition.

Definition 1. Let $\beta$ be ${k}$-dimensional vector and $\sigma > 0$, if the joint distribution of $(\beta, \sigma)$ has the following form

where C is a proportionality constant, then we have that $(\beta, \sigma)$ follows the k-dimensional distribution of $NIG_k (\mu, \Lambda, a, b)$, that is, $f\left(\sigma \right)$$\sim$$IG\left(a, b \right)$, $f\left(\left. \beta \right|\sigma \right)\sim {{N}_{k}}\left(\mu, \sigma {{\Lambda }^{-1}} \right)$.

Let $\hat{\beta}_{\tau} = {{({{X}^{*T}}{{\Sigma }^{-1}}{{X}^{*}})}^{-1}}{{X}^{*}}{{\Sigma }^{-1}}{{Y}_{\tau}^{*}}$, the likelihood (2.8) can be reformulated as:

where ${{\mu }_{\tau }} = {{\hat{\beta }}_{\tau }}$, $\Lambda = {{X}^{*T}}{{\Sigma }^{-1}}{{X}^{*}}$, $a = \frac{n-k-2}{2}$ and $b_{\tau} = \frac{1}{2}{{[{{Y}^{*}}-{{X}^{*}}{{{\hat{\beta }}}_{\tau }}]}^{T}}{{\Sigma }^{-1}}[{{Y}^{*}}-{{X}^{*}}{{{\hat{\beta }}}_{\tau }}]$.

According to Definition 1, the likelihood function (2.8) can be represented as the $NIG({{\mu }_{\tau }}, \Lambda, a, b_{\tau})$ distribution,

Now, for the sake of computational convenience, a conjugate NIG prior distribution is needed to construct our Bayesian model of streaming data. According to the study of Chu et al. ^[31], a conjugate prior for $(\beta, \sigma)$ with a modification of Zellner's informative $\it g$-prior in quantile regression could be expressed as:

Based on the above prior distributions of $\beta$ and $\sigma$, a joint $g$-NIG prior distribution can be obtained by

where ${{\mu }_{0}} = {{0}_{k}}$, ${{\Lambda }_{{{g}_{0}}}} = \frac{{{X}^{*T}}{{\Sigma }^{-1}}{{X}^{*}}}{g}$, ${{a}_{0}} = 0$, ${{b}_{0}} = 0$, $g > 0$ is a known parameter and can be set freely.

By (2.11) and (2.14), we have the following hierarchical model.

2.3.   The full conditional distributions and the Gibbs sampler

Given the informative $g$-prior distribution, the model becomes computationally easier to solve since the full conditional distributions now are

(1) Sample $v_i$ from

where, $\bar{\xi}_i = \frac{(y_i-x_i^T\beta)^2+(\beta^T x_i x_i^T \beta)/g}{2\sigma}$, $\bar{\zeta} _i = \frac{1}{2\sigma}, i = 1, 2, \cdots, n.$

(2) Sample $\beta$ from

where, ${{\bar{\mu }}_{\tau }} = {{[(1+\frac{1}{g}){{X}^{*T}}{{\Sigma }^{-1}}{{X}^{*}}]}^{-1}}{{X}^{*T}}{{\Sigma }^{-1}}Y_{\tau }^{*}$, $\bar{\Lambda } = (1+\frac{1}{g}){{X}^{*T}}{{\Sigma }^{-1}}{{X}^{*}}$.

(3) Sample $\sigma$ from

where, $a = \frac{3n+k}{2}$, $b = \frac{1}{2}[{{(Y_{\tau }^{*}-{{X}^{*}}\beta)}^{T}}{{\Sigma }^{-1}}(Y_{\tau }^{*}-{{X}^{*}\beta })+{{\beta }^{T}}\frac{{{X}^{*T}}{{\Sigma }^{-1}}{{X}^{*}}}{g}\beta +2\tau (1-\tau)\sum\limits_{i = 1}^{n}{{{v}_{i}}}].$

Here, $GIG(0, {{{\bar{\xi }}}_{i}}, {{{\bar{\zeta }}}_{i}})$ is the generalized inverse Gaussian distribution with probability density function $x^{-1}exp\{\frac{1}{2}(\bar{\zeta}x^{-1}+\bar{\xi}x)\}$ for $x > 0$.

For the convenience of calculation and expression, based on the posterior distribution of $\beta$ and $\sigma$, we have a joint posterior distribution of $\pi(\left. \beta, \sigma \right|Y_{\tau }^{*}, {{X}^{*}}, V)$ as follows:

where, ${{\bar{\mu }}_{\tau }} = {{[(1+\frac{1}{g}){{X}^{*T}}{{\Sigma }^{-1}}{{X}^{*}}]}^{-1}}{{X}^{*T}}{{\Sigma }^{-1}}Y_{\tau }^{*}$, $\bar{\Lambda } = (1+\frac{1}{g}){{X}^{*T}}{{\Sigma }^{-1}}{{X}^{*}}$, $\bar{a} = \frac{3}{2}n$, ${{\overline{b}}_{\tau }} = \frac{1}{2}{Y_{\tau }^{*}}^T{{\Sigma }^{-1}}Y_{\tau }^{*}-\frac{1}{2}\bar{\mu }_{\tau }^{T}\bar{\Lambda }{{\bar{\mu }}_{\tau }}+\tau (1-\tau)\sum\limits_{i = 1}^{n}{{{v}_{i}}}$.

3.   Bayesian quantile regression for streaming data

Different from traditional regression datasets, our focus in this paper is on streaming data sets. Assume $\{D_m, m = 1, \cdots, M, \cdots\}$ are aggregated streaming data, $D_m = \{(y_{mi}, \mathit{\boldsymbol{x}}_{mi}), i = 1, \cdots, n_m\}$ are the streaming data of $\it m$-th batch with sample size $n_m$, and the total sample size up to $M$ is $N_M = \sum_{i = 1}^{M} n_i$. To construct the posterior distribution for streaming data, the NIG multiplier operator defined in Chu et al. ^[31] will be used, as presented in the following proposition.

Proposition 1. The multiplication of H $k$-dimensional distributions $NIG_k(\mu_h, \Lambda_h, a_h, b_h)$ is also a NIG distribution, that is

where $\mu = (\sum_{h = 1}^H \Lambda_h)^{-1} (\sum_{h = 1}^H \Lambda_h \mu_h)$, $\Lambda = \sum_{h = 1}^H \Lambda_h$, $a = \sum_{h = 1}^H a_h+\frac{(H-1)(k+2)}{2}, b = \sum_{h = 1}^H b_h+\frac{1}{2} \sum_{h = 1}^H(\mu_h-\mu)^{T} \Lambda_h(\mu_n-\mu)$.

We begin with the $\it m$-th batch of streaming data $D_m$. Let ${Y}_{\tau, m}^{*}$ and ${X}_{ m}^{*}$ be the data in $\it m$-th batch $D_m$, ${\Sigma }_m$ be an $n_m \times n_m$ diagonal matrix, $\hat{\beta}_{\tau, m} = {{\left({{X}_m^{*T}}{{\Sigma }_m^{-1}}{{X}_m^{*T}} \right)}^{-1}}{{X}_m^{*}}{{\Sigma }_m^{-1}}{{Y}_{\tau, m}^{*}}$. From (2.10), the likelihood function can be derived for the $\it m$-th batch of streaming data as follows:

where, ${\mu }_{\tau, m} = \hat{\beta}_{\tau, m} = {{({{X}_m^{*T}}{{\Sigma }_m^{-1}}{{X}_m^{*}})}^{-1}}{{X}_m^{*}}\Sigma_{m}^{-1}{{Y}_{\tau, m}^{*}}$, $\Lambda_m = {{X}_m^{*T}}{{\Sigma }_m^{-1}}{{X}_m^{*}}$, $a_m = \frac{n_m-k-2}{2}$, $b_{\tau, m} = \frac{1}{2}{{[{{Y}_{\tau, m}^{*}}-{{X}_m^{*T}}{{{\hat{\beta }}}_{\tau, m }}]}^{T}}{{\Sigma }_m^{-1}}[{{Y}_{\tau, m}^{*}}-{{X}_m^{*T}}{{{\hat{\beta }}}_{\tau, m }}]$.

Subsequently, the posterior distribution of the streaming data can be obtained by multiplying the NIG prior of the stream data with the likelihood function of the $m$-th batch data. To implement our BSQR method, the prior distribution of the $m$-th batch of streaming data is introduced in the following definition.

Definition 2. The prior distribution of the $\it m$-th batch is defined by $NIG(\mu_{p, \tau, m}, \Lambda_{p, m}, a_{p, m}, b_{p, \tau, m})$, which consists of two parts, the adjusted posterior distribution of historical data and the $\it g$-NIG prior of current data, that is

where, $NIG^{(m-1)}_k$ denotes the posterior distribution of the $(m-1)$-th batch for streaming data,

Remark 1. For streaming data $D_1, D_2, \cdots D_m, \cdots$, if $n_1 \ge \frac{k+2}{3}$, where $n_1$ is the sample size of the first batch of streaming data, the prior distribution $NIG(\mu_{p, \tau, m}, \Lambda_{p, m}, a_{p, m}, b_{p, \tau, m})$ exists.

Through the prior distribution of stream data $NIG(\mu_{p, \tau, m}, \Lambda_{p, m}, a_{p, m}, b_{p, \tau, m})$, we can get the posterior distribution of the streaming data. Theorem 1 describes the renewability of the posterior distribution of streaming data obtained through the NIG multiplier operator.

Theorem 1. Consider a linear QR model with streaming data observations $D_1, D_2, \cdots D_m, \cdots$. Setting the prior distribution of the initial streaming data as the g-NIG prior $NI{{G}_{k}}\left({\mu }_{0}, {{\Lambda }_{{{g}_{0}}}}, {a}_{0}, {b}_{0} \right)$. When the $m$-th batch of streaming data is received, the posterior distribution of the streaming data $N I G_k^{(m)}(\bar{\mu}_{\tau, m}, \bar{\Lambda}_m, \bar{a}_m, \bar{b}_{\tau, m})$ can be expressed as

where, $N I G_k^{(L)}(\mu_{\tau, m}, \Lambda_m, a_m, b_{\tau, m})$ is the likelihood function for the $m$-th block of data nodes,

$\bar{\Lambda }_m = \bar{\Lambda}_{m-1}+\Lambda_{m, 0}+\Lambda_m$ $= (1+\frac{1}{g})\sum_{i = 1}^{m}{{X}_i^{*T}}{{\Sigma }_i^{-1}}{{X}_i^{*}}$,

$\bar{a}_m = \bar{a}_{m-1}+a_0+a_m+\frac{k+2}{2} = \frac{3}{2}n_m+a_0$,

${{\overline{b}}_{\tau, m }} = \bar{b}_{\tau, m-1}+ b_0 +b_{\tau, m} + \frac{1}{2}(\bar{\mu}_{\tau, m-1}-\bar{\mu}_{\tau, m})^T\bar{\Lambda}_{m-1}(\bar{\mu}_{\tau, m-1}-\bar{\mu}_{\tau, m}) +\frac{1}{2}\bar{\mu}^T_{\tau, m}\Lambda_{m, 0} \bar{\mu}_{\tau, m}$+$\frac{1}{2}(\mu_{\tau, m}-\bar{\mu}_{\tau, m})^T\Lambda_m(\mu_{\tau, m}-\bar{\mu}_{\tau, m})$ $= \frac{1}{2}\sum_{i = 1}^mY_{\tau, i }^{*T}\Sigma_{i}Y_{\tau, i }^{*}-\frac{1}{2}\bar{\mu }_{\tau, m }^{T}\bar{\Lambda }_m{{\bar{\mu }}_{\tau, m }}+\tau \left(1-\tau \right)\sum\limits_{i = 1}^{n}{{{v}_{i}}}+b_0$.

According to Theorem 1, the update of parameters in the posterior distribution of the streaming data only requires the parameters of the posterior distribution of historical data. Specifically, after receiving the streaming data, we only need to calculate and store the parameters of the posterior distribution $\bar{\mu}_{\tau, m}$, $\bar{\Lambda}_{m}$, $\bar{a}_m$ and $\bar{b}_{\tau, m}$. In addition, when given all the values of $v_i$, we can see from Theorem 1 that the posterior distribution parameters of streaming data are equal to the posterior distribution parameters of the entire data. The proof of Theorem 1 is provided in the appendix.

Under Theorem 1 and Proposition 1, an effective algorithm for the Bayesian streaming quantile regression (BSQR) is developed in Algorithm 1.

4.   Simulation study

In this section, the proposed BSQR will be illustrated by simulated data.

We consider datasets in different cases and generate the streaming datasets $D_1, D_2, \cdots, D_M$ up to batch $M$, where $D_m = \{(y_{mi}, \mathit{\boldsymbol{x}}_{mi}):1\le i\le n_m \}$.

where, the elements of $\mathit{\boldsymbol{x}}_{mi}$ are independently and identically distributed as the standard normal distribution, and the vector $\mathit{\boldsymbol{\beta}}_0 = (1, 3, 2, 10, 4, 3, -1, 4, 5, 0)^T$. For the random error $e_{mi}$, in order to illustrate the robustness, the following three cases are considered.

Case 1. $e_{mi}$ follows the asymmetric Laplace distribution.

Case 2. $e_{mi}$ follows the standard normal distribution.

Case 3. $e_{mi}$ follows $t(3)$ distribution.

For a comprehensive illustration, we compare our method with four competing estimators in the current literature.

(1) The oracle Bayesian quantile regression (BQR) obtained by processing the entire data once.

(2) The renewable quantile regression (RQR) proposed in Wang et al. ^[30].

(3) The renewable SQR (RSQR) estimator proposed in Jiang et al. ^[37].

(4) The online linear estimator for the QR (OLEQR) proposed in Chen et al. ^[28].

The first competing estimator BQR is used to show that our method is statistically equivalent to the oracle one and verify the theoretical results of Theorem 1 for the situation of finite sample. The second to fourth estimator RQR is used to compare our method with the method of streaming data in the frequentist framework, which makes our simulation experiments more comprehensive.

For the selection of hyperparameter g in g-NIG prior, Smith and Kohn^[33] found that choosing $g$ between 10 and 1000 can achieve good results. Meanwhile, they set $g = 100$ as the final value. Dao and Wang ^[34] also used $g$-NIG priors and found that their numerical results for different choices of $g$ within the range of $100 \le g \le 1000$ yielded very similar results. For our model, we also validated the selection of g using simulation experiments, and the results are shown in Figure 1. It can be observed that the choice of g has little impact on MSE. Therefore, we follow the suggestions of Smith and Kohn and set $g = 100$, a choice adopted by various researchers, such as Lee et al.^[35] and Chen et al.^[36].

To examinate the perfomance of the BSQR method, we fix $n$ with varying batch size $n_m$. Specifically, $M$ streaming datasets with $n = 10^4$ observations are generated independently, where batch number $M$ is between 50 and 100, with each batch having $\frac{10000}{M}$ observations. We repeat the experiment 500 times and define $\hat{\beta}_{j(s)}$ as the estimator of $\beta_j$ in the $\it s$-th replication. In order to compare estimation accuracy, the mean square errors (MSE) $MSE(\hat{\beta}_j) = \frac{1}{500}\mathop{\sum}\limits_{s = 1}^{500}(\hat{\mathit{\boldsymbol{\beta}}}_{j(s)}-\mathit{\boldsymbol{\beta}}_{0j})^2, j = 1, 2, \cdots, 10$, are calculated for the three methods. The results of case 1 are shown in Tables 1 and 2.

In practice, the error term ${e}_{mi}$ may not be distributed as AL distribution. Therefore, Cases 2 and 3 are also simulated when the error term follows t-distribution with a degree of freedom 3, and standard normal distribution, respectively. The results are shown in Tables 3–6.

From the simulation results, we have several findings. First, our streaming data Bayesian quantile regression method can always estimate the parameters consistently under three error distributions. In the meantime, the estimation results are not affected by the data batch size, since the MSE values of our method are almost in the same scale as the batch number increases from 50 to 100. Second, under ALD error distribution, our method performs better than the other three methods, and also performs comparably with the oracle Bayesian quantile regression, as the BQR and BSQR methods are very close. Third, under other error distributions, our method performs worse than RQR and RSQR methods, but still performs better than the OLEQR method.

Additionally, we consider an interesting scenario 2 where streaming datasets arrive at a high speed and the batch size is small. For convenience, we fix the batch size $n_1 = \cdots = n_m = 100$, and let n increases from $10^5$ to $10^6$. The results are shown in Table 7 with $\it ALD$ error distribution. The results for other quantile levels and error distributions are similar and thus omitted. We can see from Table 7 that the performance of our method improves (the MSE values all decrease) as the number of batches increase from 1000 to 10000, and it always works comparably with the oracle Bayesian quantile regression estimator. These results imply that our method is not affected as more data streams are processed; in other words, it is robust to the number of batches.

For Algorithm 1, we conduct an analysis of its computational complexity. By analyzing the computational process of Algorithm 1, we determine that its computational complexity is $O((n_mk)^2)$. Furthermore, we record its running time and compare it with the distributed algorithm. For the convenience of recording and computational efficiency, we fix the batch size $n_1 = \cdots = n_m = 100$, and let $m$ increases from $10$ to $100$. All results are based on 1000 draws obtained from the Gibbs samplers. Through this comparative analysis of running times, we observe that the streaming data algorithm has better computational efficiency. The comparison result is shown in Figure 1, where the horizontal axis represents the number of batches and the vertical axis represents the running time unit in minutes. Additionally, by comparing Tables 1–7 and Figure 2, it can be observed that our method is similar to the full data method, but our method has better effiency than the distributed algorithm. Therefore, the BSQR method can be used as an alternative to distributed algorithms to improve computational efficiency, especially when the number of data batches is large.

5.   Real data analysis

Example 1. The online news popularity dataset. In this example, we illustrate our Bayesian streaming quantile regression method by using the Online News Popularity dataset. The dataset is from the UCI machine learning repository. This is a large-scale dataset with $N = 39,640$ observations. In our analysis, the number of shares in social networks is set as the response variable, and the perspective of web environment factors such as positive words, topic keywords, etc., are chosen as 7 predictors. For comprehensiveness, five different quantile levels from small to large are considered, i.e., $\tau = 0.1, 0.3, 0.5, 0.7$, and 0.9.

We fit our method to the above dataset by implementing Algorithm 1. In each $\tau$, we partition the dataset into 100 subsets with the size of $n_m = 396$ for $m = 1, \cdots, 99$ and $n_{100} = 436$. All results are based on 15,000 draws obtained from the Gibbs samplers with a burn-in of 5000 iterations. The estimated coefficients and posterior standard deviations at the specified quantile levels are presented in Table 8.

We observe that the number of positive words is negatively correlated with the popularity, while length, num-img, num-video, and num-negitive have positive impacts across all quantiles. Additionally, as the quantile increases, the impact of title sentiment polarity first increases and then decreases, which means that the title of a popular news cannot be too polarized. This empirical study shows that the BSQR method we proposed is helpful in investigating the impact of different factors on news popularity in streaming data scenarios.

Example 2. The Beijing multi-site air-quality dataset. This example illustrates the streaming data algorithm for BSQR by using the Beijing multi-site air-quality dataset. The data were collected from the Beijing Municipal Environmental Monitoring Center, which is a reputable institution responsible for monitoring air quality, meteorological parameters, and other environmental factors in Beijing. The dataset includes measurements of pollutants such as PM$_{2.5}$, PM$_{10}$, NO$_2$, SO$_2$, CO, and O$_3$, as well as meteorological variables like temperature, humidity, wind speed, and wind direction. The data were recorded at multiple monitoring stations across the city, providing a spatially diverse dataset. The dataset is publicly available at https://www.bjmemc.com.cn/, spanning the period from July 1, 2017 to March 30, 2023. The dataset includes 430,524 hourly air pollutant data points from 12 district-controlled air quality monitoring stations in Beijing. The reasons for using air quality data are as follows: (1) Air quality data necessitates frequent real-time updates to accurately depict the current state of air quality. Streaming data frameworks are capable of rapidly processing and providing feedback on real-time data. (2) Air quality monitoring stations continuously and extensively generate data. These stations are dispersed across various locations, with each site persistently producing a substantial volume of data. Streaming data frameworks can efficiently handle this continuous data flow.

In this study, we explore the relationship between the PM$_{2.5}$ concentration (ug/m$^3$) and ten variables in Table 9. For model diagnostics, we use the AIC method for model selection. Specifically, we collected the historical Beijing multi-site air-quality dataset from March 2013 to February 2017, and calculated the AIC values under three different models based on the historical data, namely the classic linear regression model and the quantile regression model with $\tau = 0.1$ and $0.9$. The AIC results are shown in Figure 3. From Figure 3, we can observe that the model with 10 variables is the best model.

Because this data was collected by air quality monitoring sites in chronological order, we set the number of batches $b = 48$ based on the number of months. This air quality dataset and model have been analyzed in Jiang and Yu ^[37], which also focused on quantile regression for streaming data.

We fit the streaming data to the Bayesian quantile regression model by implementing Algorithm 1 for five quantile levels, specifically $\tau =$ 0.10, 0.30, 0.50, 0.70, and 0.90. In each scenario, the entire set of observations is divided into 48 subsets. Subsequently, the informative g-prior is assigned with $g$ = 100. All reported results are based on 15,000 draws obtained from the Gibbs samplers, with a burn-in period of 5000 iterations. Table 10 and Figure 4 present the estimated coefficients for the specified quantile levels.

At different quantiles, the coefficients of five air factors are all significant and positive, indicating that these air pollutants have a positive impact on PM$_{2.5}$ concentration. As the quantile $\tau$ increases, the estimated coefficients of SO$_2$, NO$_2$, and O$_3$ decrease, which means that the impact of SO$_2$, NO$_2$, and O$_3$ on PM$_{2.5}$ concentration gradually decreases with the increase of PM$_{2.5}$ concentration. Additionally, the estimated coefficients of PM0 and CO increase with increasing quantile $\tau$, indicating that as the concentration of PM$_{2.5}$ increases, the impact of PM$_{10}$ and CO is increasing.

For the other five meteorological factors, TEMP, PRES, DEWP, RAIN, and WSPM, RAIN and WSPM have significant and negative coefficients at each quantile, indicating that they have a negative impact on PM$_{2.5}$ concentration, and this effect become weaker as PM$_{2.5}$ concentration increases. The coefficient of PRES is not significant at the 0.9 quantile, indicating that air pressure has no significant effect on PM$_{2.5}$ at high concentrations. The TEMP coefficient is significantly negative at five quantiles, and the absolute value of the coefficient gradually increases, meaning that as the concentration of PM$_{2.5}$ increases, the impact of temperature becomes stronger. The coefficients of DEWP at each quantile are significantly positive, and as the quantile increases, the coefficients show an increasing trend. This indicates that as the concentration of PM$_{2.5}$ increases, the impact of DWEP increases.

Compared with the method in ^[37], we added three new covariates, PM$_{10}$, O$_3$, and RAIN, making the analysis of the real data more comprehensive. In addition, the conclusion of BSQR method indicates that the estimated coefficients of SO$_2$ and NO$_2$ decrease as quantity $\tau$ increases, and the estimated coefficients of WSPM are negative. These findings are different from the conclusions of ^[34], but are consistent with the research results of PM$_{2.5}$ in Li^[38] and Zhang and Zhang^[39].

6.   Conclusions

This paper extends the traditional frequency QR. This extension involves the utilization of an ALD-based likelihood function, a conjugate NIG prior, and an adjusted NIG posterior. The streaming data method we proposed does not require storing historical data; it only needs to save the parameters of the posterior distribution and the current data. Furthermore, theoretical results establish the equivalency between the proposed posterior distribution for streaming data and the one computed using the full data. Both simulation studies and real data evaluations affirm the strong performance of the new approach proposed in this paper.

Real data analysis of Example 1 shows that the popularity of news is related to factors such as the number of positive words, images, video content, keywords and length. Among them, the number of positive words has a negative impact on the popularity of news. In Example 2, air pollutants and meteorological factors have varying impacts on PM$_{2.5}$ concentration, with the influence of some factors changing as PM$_{2.5}$ concentration increase. Specifically, the impact of air pollutants SO$_2$, NO$_2$, and O$_3$ decreases, while that of PM$_{10}$ and CO increases. Among meteorological factors, the influence of RAIN and WSPM decreases, and the impact of TEMP and DEWP becomes more pronounced at high PM$_{2.5}$ concentrations.

This study reveals the degree and trend of the impact of different air pollutants and meteorological factors on PM$_{2.5}$ concentrations, which helps to more accurately understand the complex mechanisms of PM$_{2.5}$ formation and provides a scientific basis for air quality management. The following air quality management strategies are recommended: (1) Prioritize the control of SO$_2$, NO$_2$, and O$_3$ emissions, especially in areas with lower PM$_{2.5}$ concentrations, as these pollutants have a significant impact on PM$_{2.5}$; (2) in areas with higher PM$_{2.5}$ concentrations, strengthen the control of PM$_{10}$ and CO to reduce their contribution to PM$_{2.5}$ concentrations; (3) strengthen the study of air pressure to determine its impact at different PM$_{2.5}$ concentrations, providing a more comprehensive basis for management strategies.

Nonetheless, there are several issues requiring future research. First, in conjunction with penalty functions such as LASSO ^[40] and SCAD^[41], it is necessary to further explore the high-dimensional regularization for streaming data. Specifically, there are several versions of regularization that could be used^[42]. In the frequency framework, Ma et al.^[43] proposed an online updating coordinate descent algorithm, and a tuning lasso parameter selection was also suggested. Wang et al. ^[45] proposed a provable online feature selection algorithm that utilizes the online leverage score. For the Bayesian framework, we will explore Bayesian variable selection strategies in the streaming data context in our subsequent work. Based on current research, we aim to propose more accurate and efficient regularization methods for streaming data. Second, the new approach focuses on Bayesian linear QR models, and the issue of concept drift ^[44,46] is not considered. Thus, one can further consider more complex Bayesian streaming data models like Bayesian deep neural networks, and invest the concept drift problem in the online BQR setting. Lastly, the new algorithm is geared to independent data, while in practical scenarios, dependent or non-stationary data ^[47] are more commonly used, including time series data or longitudinal data. Hence, it will be interesting to develop the corresponding algorithms for such kinds of datasets.

Acknowledgments

This research was funded by Natural Science Basic Research Program of Shaanxi (Program No. 2023-JC-QN-0059).

We would like to express our sincere gratitude to the editor and the reviewers for their valuable time, efforts, and expertise in reviewing our manuscript.

Conflict of interest

The authors declare no conflict of interest.