Demand Forecasting with Bayesian Methods

Online Regression

In a typical Bayesian regression problem, we have a setting with independent and identically distributed training data with future observations drawn from the same population. In this case we can simply place a prior over the weights and variance, evaluate the data likelihood and find the posterior through exact or approximate inference to determine a stable model from which we can make predictions. Here we are presented with a time series of features and inventory data and want our model to flexibly model trends that depend on our time step. Online regression is a naive baseline for this problem in which at time step t, we take the posterior from time step t − 1 as our prior and update the model using only the likelihood over the most recent batch of data. Define our original priors as \(p(\theta_0)\), then our model updates occur sequentially:

\[\begin{align*} p(\theta_{t}) = p(\theta_{t} \mid y_{t-1}, x_{t-1},\theta_{t-1}) &\propto p( y_{t-1}\mid x_{t-1},\theta_{t-1})p(\theta_{t-1})\\ p(\theta_{t+1} \mid y_{t}, x_{t},\theta_{t}) &\propto p( y_{t}\mid x_{t},\theta_{t})p(\theta_{t}) \end{align*}\]

While a simple regression may seem underpowered for this task, it provides an excellent baseline against which we can compare more complex approaches. The conjugacy of the Gaussian model makes it extremely efficient as we add new data points with easy interpretation. However, its failure to model dependency between data points may introduce bias that could be avoided with more complex approaches.

The model as implemented closely follows the approach outlined in Christopher Bishop’s Pattern Recognition and Machine Learning. Given model \(y_i = w_ix_i^T + \varepsilon_i\), we give likelihood and prior distributions as \(p\left(y_{i} \mid x_{i}, w_{i}\right)=\mathcal{N}\left(w_{i} x_{i}^{\top}, \beta^{-1}\right)\) and \(p\left(w_{0}\right)=\mathcal{N}\left(m_{0}, S_{0}\right)\) respectively, where \(\beta\) is a precision hyperparameter, \(m_0 = \mathbf{0}\) and \(S_0 = diag(\alpha^{-1})\) for hyperparameter \(\alpha\). Our posterior is analytically derived as \[\begin{align*} p\left(w_{i+1} \mid w_{i}, x_{i}, y_{i}\right)=\mathcal{N}\left(m_{i+1}, S_{i+1}\right)\\ S_{i+1}=\left(S_{i}^{-1}+\beta x_{i}^{\top} x_{i}\right)^{-1}\\ m_{i+1}=S_{i+1}\left(S_{i}^{-1} m_{i}+\beta x_{i} y_{i}\right) \end{align*}\]

For model criticism, Bishop recommends evaluating the marginal likelihood, or model evidence, which he derives analytically. For M features and N observations, the log marginal likelihood can be expressed as \[\begin{align*} \log p({y} \mid \alpha, \beta)=\frac{M}{2} \log \alpha+\frac{N}{2} \log \beta-E\left({m}_{N}\right)-\frac{1}{2} \log \left|{S}_{N}^{-1}\right|-\frac{N}{2} \log 2 \pi \end{align*}\] Therefore we can directly use this expression to select for optimal hyperparameters \(\alpha,\beta\) to improve the model fit.

Gaussian Processes Regression

I look to Gaussian processes (GP) regression to correct for the deficiencies of the regression model. While the prior model considered datapoints to be independently drawn from a stationary distribution, Gaussian processes explicitly model the covariance relations between time steps. Here we consider our feature vectors \(x \in R^d\) over T time steps \(x_1,...,x_T\) and each observation of inventory used \(y_i = f(x_i)\). If we assume the functions \(f\) to be given by a prior Gaussian process \(f \sim GP(\mu, K)\) for some mean \(\mu\) and covariance matrix \(K\) and \(y_i \sim N(f(x_i), \sigma^2)\) then we obtain posterior \(p\left({f} \mid x_{n}, y_{n}\right) \propto {N}\left({f} \mid {\mu}^{\prime}, {K}^{\prime}\right)\) for \({K}^{\prime}=\left({K}^{-1}+\sigma^{-2} {I}\right)^{-1}, \; {\mu}^{\prime}={K}^{\prime}\left({K}^{-1} {\mu}+\sigma^{-2} {y}\right)\). Similarly, we can obtain an exact posterior predictive distribution given by \[\begin{align*} y_{T+1} | x_{T+1}, x_i, y_i &\sim N(m_{T+1}, v_{T+1}) \\ m_{N+1} &=\mu\left({x}_{N+1}\right)+{k}^{\top}\left({K}+\sigma^{2} {I}\right)^{-1}({y}-{\mu}) \\ v_{N+1} &=K\left({x}_{N+1}, {x}_{N+1}\right)-{k}^{\top}\left({K}+\sigma^{2} {I}\right)^{-1} {k}+\sigma^{2} \end{align*}\]

GP regression offers additional modeling decisions that improve its attractiveness for our use case. First, there are a variety of potential kernels for expressing the relationship between feature vectors over time. Our choice determines how closely we model regular periodicity and linear trends in the data.

Additionally, since GP regression is closely related to time series methods like ARIMA, we can eschew using features at all and treat the blood inventory as a standalone data series. Why might we prefer this approach? The other regression approaches obtain a posterior predictive distribution by conditioning on \(x_{N+i}\) for \(i \in [1,7]\), feature vectors which are supposed to contain data like patient counts and other measures that aren’t actually available for future days. While features could be extrapolated from current data, a GP regression approach can condition simply on the time step and avoid some of these messier modeling decisions.

Online Regression

The regressions were tested under two scenarios: for each time period, we could either fit the model using data since inception or look back over a fixed window. The data since inception approach was too constraining and performed slightly worse in evaluations, leading the fixed window approach to be used for comparisons. The model was fit over 7 day increments and produced predictions for 7 days into the future.

We can understand the performance of the model by looking at the results more closely. First, while the online regression approach does not explicitly model a time series, one-hot encoded day of the week features are included as input to the model. Therefore it still can capture some of the cyclicality of blood platelet demand.

The online regression was fit over a grid of potential hyperparameter values for \(\alpha, \beta\), with the model with the highest marginal likelihood selected for final predictions. It was sensitive to extreme values of each hyperparameter, but otherwise relatively robust to misspecification within a neighborhood of values.

Gaussian Processes Regression

GP regression offers a variety of modeling decisions, but perhaps the most important is the choice of kernel. This decision directly impacts the model’s ability to reflect the true relationship between observations. In “Automatic model construction with Gaussian processes” David Duvenaud outlines a decision process for choosing the best kernel given a dataset. I experimented with permutations of additive kernels, combining Exponential Quadratic, Linear, Periodic, and Locally Periodic varieties. The Locally Periodic kernel was of particular interest as the data shows consistent day of the week seasonality but with magnitudes that vary over time.

The model was fit in TensorFlow Probability by maximizing the marginal likelihood of the model parameters and the observation noise variance over the training data via gradient descent. The GP regression model was then fit on the data using the optimized kernel before posterior samples were drawn. Over all kernel experiments, a combination kernel of Exponential Quadratic and Locally Periodic varieties produces the best results.

Left: Mean and credible intervals drawn from the posterior predictive distribution. Right: sample paths drawn from the posterior predictive distribution.

The model was then fit for periods increasing in 7 day increments from the data series inception, with predictions made for the following 7 day increment. For example at day 100, data from days 0-100 are taken as the training set, with the mean and standard deviation from the posterior predictive distribution forming confidence intervals for the prediction period. Alternatively, samples could be drawn from the posterior predictive distribution to simulate different paths created by the Gaussian process.

Left: Samples drawn from an Exponential Quadratic / Locally Periodic kernel with parameters learned from fewer than 100 datapoints. Middle: The same kernel with parameters learned from all datapoints. Right: First 100 observations of platelets used.

We can make further sense of GP regression’s underperformance by looking at the kernels learned themselves. Comparing kernels fit in the training process along with the true observations, we can see that given little data, the kernel is unable to adapt to the variation in the data. A kernel with parameters trained on the full dataset shows much more reasonable variation, with both global and local trends. However, comparing this kernel to the observed variation, we can see it still has a long way to go to match the ground truth. Additional data would go a long way to improving this method, but there may also be more expressive kernels than those explored that could drastically improve performance.

While the hyperparameter search was partially limited by compute time, GP regression’s results indicate the method could show promise in this application. While it achieved higher error than other methods, its performance could improve if we included additional features, some of which might be known ahead of time by hospital administrators. Additional variance control, such as full marginalization of the model’s hyperparameters, would also increase its effectiveness.