Using infinite server queues with partial information for occupancy prediction

Abstract

Motivated by demand prediction for the custodial prison population in England and Wales, this paper describes an approach to the study of service systems using infinite server queues, where the system has non-empty initial state and the elapsed time of individuals initially present is not known. By separating the population into initial content and new arrivals, we can apply several techniques either separately or jointly to those sub-populations, to enable both short-term queue length predictions and longer-term considerations such as managing congestion and analysing the impact of potential interventions. The focus in the paper is the transient behaviour of the $M_t/G/\infty$ queue with a non-homogeneous Poisson arrival process and our analysis considers various possible simplifications, including approximation. We illustrate the approach in that domain using publicly available data in a Bayesian framework to perform model inference.

Keywords:

• Queuing
• decision support systems
• government
• health services

1. Introduction

This work is motivated by the problem of predicting short and longer-term implications of policy changes on the custodial elements of the prison system in England and Wales. The model described here was developed following consultation with the Ministry of Justice (MoJ) to add to their methods of forecasting the prison population, to help analyse the implications of changing external factors accounting for the prison population such as government guidelines and sentencing policies, and to consider the uncertainty involved in the model and its predictions.

The nature of the prison system is such that arrivals can’t be turned away, hence infinite server queueing models are directly applicable as they support the assumption of no queueing delay for service. While such models have been widely used in modelling service systems, including call centres (Ibrahim et al., 2016), hospital staffing (Pender, 2016), software reliability (Huang & Huang, 2008) and insurance claims (Cheung et al., 2019), the assumptions relevant to our scenario lead us to consider some less well-known results from the queueing systems literature and discuss how they can be useful in this setting.

In this domain, matters of capacity management and overcrowding at individual institutions have to be handled by adjustments involving medium term system-wide considerations such as sentencing patterns, parole guidelines and the use of community service arrangements. This intent to support policy makers leads to a focus on the strategic level of decision support (Grieco et al., 2021; Hulshof et al., 2012), and to some considerations on the development of models expressed in terms that are interpretable by policy makers and can enable “what if?” studies, including quantification of the impact of prospective policy change (Bravo et al., 2019; Dong & Perry, 2020; Kegel et al., 2017; Petris et al., 2009; Tuominen et al., 2022).

The value of decision support tools that can analyse the impact of interventions linked to policy change has been demonstrated in care pathways in health services (Demir & Southern, 2017). While modelling flows through a prison system has not been widely studied, there has been some work using queueing models to study the relationship between policy changes and prison occupancy. For example, Usta and Wein (2015) used a queueing network model of the jail system and a simulation approach to study the optimal mix of pretrial release and forms of sentencing to minimise the amount of recidivism, subject to constraints on the available prison occupancy; and Master et al. (2018) used a $M/M/c/c$ queueing system to assess performance of alternative pretrial release policies, and of sentences with a split of custodial and supervised outcomes.

Using public data, we study a single phase of a prisoner’s journey through the prison system, as it is well out of scope to model the full system. Our analysis is both informed by, and limited by, the availability of service system data: the size of the prison population is collected and recorded on a monthly basis, but the time served and the remaining length of stay (service times) of those individuals are not. Indeed for custodial sentence admissions, even within offence types (the 12 main offence groups are listed in Table 1), sentence lengths differ, and there is a difference between a court imposed sentence length and the actual service time, so the model cannot assume remaining service times for those individuals present at a given time, even if their formal sentence lengths are known. This contrasts with, for example, work on bed demand in an intensive care unit that considered existing patients as well as arrivals, but could use knowledge of the length of stay of existing patients (Pagel et al., 2017).

Queueing model analyses typically assume the system begins in an empty state (Li et al., 2019) and then, assuming a Poisson arrival process, the well-known results of Eick et al. (1993) mean that the queue length exhibits a Poisson distribution with a mean derivable from properties of the service time distribution. The scenario motivating the work in this paper involves the less well-studied situation where there are n > 0 individuals initially present and where the elapsed time at t = 0 of each individual is not known (Goldberg & Whitt, 2008; Korhan Aras et al., 2017; Weber, 2005). In this case, the departure process is no longer Poisson, but the queue dynamics can be analysed as a combination involving those already in service at t = 0 and those who subsequently arrive, and the departure process is the superposition of a binomial and a Poisson process (Weber, 2005). This separation into initial content and new input allow these sub-populations to be analysed jointly or separately.

As pointed out by Pagel et al. (2017), in contrast to much research involving applications of queueing systems that rely on steady state distributions, when considering the use of models for informing short term consequences associated with operational decision making, results involving transient distributions of queueing systems become relevant.

As highlighted in a recent review of major healthcare applications of infinite server queuing models (Worthington et al., 2020), time-inhomogeneous infinite server models have been used both for predictive modelling (e.g., ward capacity planning (Bekker & de Bruin, 2010)) and for investigating the impact of policy changes (e.g., the introduction of specialised treatment centres for specified illnesses (Utley et al., 2008)). Our ambition to support the analysis of both short and longer-term policy change means we also consider time-varying arrivals.

Another key choice in formulating a queueing theory model is the assumptions regarding service time. The available data suggest a heavy tailed service time (Ministry of Justice, 2020b). Hence, we are led to the $M_t/G/\infty$ model as it is well recognised that there is an effect on the performance of the queue of the tail of the distribution as it becomes less exponential (Goldberg & Whitt, 2008).

To employ our model as a prediction tool, we must also consider the uncertainty involved as we are considering behaviour beyond the sample of the original data set used to estimate the model parameters (Gans et al., 2003). Incorporating parameter uncertainty as part of model inference is an important step to avoid overconfidence in our results (Aktekin & Soyer, 2011). We consider parameter uncertainty involved in the prediction of the system size using estimates from aggregate historical data. Techniques typically rely on summary observable data such as queue lengths, visit counts and response times, but perform inference in different ways (Spinner et al., 2015). A common approach is by Jongbloed and Koole (2001) using a Poisson mixture model and a recent survey on forecasting of the arrival process is by Ibrahim et al. (2016). We employ a Bayesian framework to perform model inference and prediction (Aktekin & Soyer, 2011; Xie et al., 2014). The Bayesian framework allows us to incorporate expert knowledge about the system into prior distributions of model parameters, which is particularly important since we have limited historical data available to us (O’Hagan, 1998).

In summary, we study the transient behaviour of the time-varying infinite server queue, $M_t/G/\infty ,$ fed by a non-homogeneous Poisson arrival process whose occupancy is observed at discrete points in time, but the time in service to that point is not known. The contributions of this paper are: (i) the novel synthesis of results from several authors about transient and stationary behaviour of the $M_t/G/\infty$ queue; and (ii) application of the approach, using Bayesian inference, to the real-world domain of prison occupancy—a domain that has not been well-studied in the literature.

The structure of this paper is as follows: in Section 2 we describe the motivating application; Section 3 outlines relevant results from the literature for analysing an observed $M_t/G/\infty$ queue; Section 4 presents a Bayesian framework for the estimation of the model parameters using domain data to illustrate the use of the model to both short-term and longer-term predictions, and includes a discussion of the mathematical assumptions of the model. In Section 5 we provide concluding remarks, with comments comparing the presented queueing theoretic approach with time-series based estimation methods. The Appendices contain further information about the prison system (Appendix A), use of the theory (Appendix B), and some examples of how analysis with our model compares with using a time-series based ARIMA model (Appendix C).

2. Motivating application: Prison occupancy

The study in this paper was produced following consultation with the Ministry of Justice (MoJ). We briefly describe the custodial elements of the prison system, that is, those that require accommodation, with more details in Appendix A. Data on the prison population is collected and managed by MoJ and the HM Prison and Probation Service (HMPPS). Statistics are regularly released as well as projections of the prison population (Ministry of Justice, 2018).

Attributes of the application domain that guided our modelling choices include: a large system of multiclass arrivals with a high offered load, no abandonments, input parameters that are subject to change and operation over a long time scale. Further, the domain data suggests the use of a non-exponential service time distribution and that a stationarity assumption is reasonable over short time frames, which allowed us to take advantage that several questions of interest have more tractable solutions under the assumptions of a stationary model.

Factors such as the conviction rate (the proportion of those arriving to court that are convicted and sentenced) and the custody rate (the proportion of those sentenced that are given custodial sentences) influence the sentenced population. Hence the size and composition of the prison population is subject to policy and legislation changes, for example, changes in Home Office (government department) resources that can affect charge rates and modifications to the sentencing guidelines (MoJ, 2019). Patterns in the published data illustrate how policy and legislative changes have had subsequent impacts on prison occupancy. Hence, from a policy maker’s perspective, being able to adjust model parameters to allow, for example, for a more serious mix of offence groups coming before the courts reflects the importance of reviewing model parameters over time.

Of course, describing the dynamics of the prison system is beyond the scope of this paper, but from a modelling perspective it requires only some reasonable assumptions to treat the flow of prisoners as an $M_t/G/\infty$ queue (Schwarz et al., 2016). Figure 1 displays our model of a simplified prisoner journey with the prison population divided into three main holding phases (displayed percentages are as of June 2019) (MoJ, 2019): (i) on remand (11%), (ii) sentenced prisoners (79%) and (iii) on recall (9%). Prisoners within the licence phase are in the community. A broader view of the prison system and its constraints would include, for example, the number of offenders on probation, staffing resources required for supervision, demands on the courts and parole hearing frequencies (Crowhurst & Harwich, 2016; Ministry of Justice, 2016, 2018).

3. An observed $M_t/G/\mathit{\infty }$ queue

Motivated by obtaining a prediction of the population given partial information, we describe results applicable to an infinite server system with Poisson arrivals and general service times observed at time $\tau >0,$ under the assumption that we have no information on when the n individuals present at this time each began their service, namely an observed $M_t/G/\infty$ queue. We define the infinite server queue in Section 3.1 and then present results for the conditional distribution for the observed queue in Section 3.2. As the n individuals initially observed at time τ complete service, this sub-population will go to zero as $t\to \infty .$ Arrivals after time τ will occur according to the Poisson dynamics described in Section 3.1. These results are used in our empirical study in Section 4.

3.1. The $M_t/G/\mathit{\infty }$ queue

The $M_t/G/\infty$ queue is a service system in which individuals arrive according to a non-homogeneous Poisson process with rate function $\lambda (t),$ for $-\infty <t<\infty$ and where the service times are independent and identically distributed (i.i.d.) and independent of the arrival process (Eick et al., 1993). Let S be the service time and denote by G its cumulative distribution function (cdf) and g its density. Assume $E[S]<\infty ,\hspace{0.17em}\lambda (t)$ is bounded and integrable and define the associated random variable S_e, the stationary excess of the service time, with cdf $G_e(t)=P(S_e\le t)$ for $t\ge 0,$ (1) $$G_e(t)=\frac{1}{E\left[S\right]}{\int }_0^t{G}^c(u)\mathit{\text{du}},$$(1) where $G^c(u)=1-G(u).$

Let Q(t) represent the number of busy servers at time t and let $m(t)=E[Q(t)].$

Theorem 3.1.

(Eick et al., 1993, Theorem 1) For an $M_t/G/\infty$ queue that was initially empty at t = – ∞, for each t, Q(t) has a Poisson distribution with mean(2) $$m(t)=E\left[{\int }_{t-S}^t\lambda (u)\mathit{\text{du}}\right]=E\left[\lambda (t-S_e)\right]E\left[S\right].$$(2)

The departure process is a Poisson process with time dependent rate function ${\lambda }^{-}(t)$, where(3) $${\lambda }^{-}(t)=E\left[\lambda (t-S)\right].$$(3)

For each t, Q(t) is independent of the departure process in the interval $(-\infty ,t].$

The impact of the service time beyond its mean can be seen from $E[S_e]=\frac{1}{2}E[S](c_s^2+1)$ where $c_s^2=\mathit{\text{Var}}(S)/E{[S]}^2$ is the squared coefficient of variation (SCV).

For $\lambda (t)=\lambda ,$ the approach to the steady state is given by (4) $$m(t)=\lambda G_e(t)E\left[S\right],$$(4) where the steady state has the insensitivity property, as $m(\infty )=\lambda E[S].$ Similarly, the transient behaviour of a stationary model that has been terminated, that is, for t < 0, $\lambda (t)=\lambda$ and zero otherwise, is $m(t)=\lambda G_e^c(t)E[S].$

In an empty system, the approach to steady state is seen from Equation (4)(4) $$m(t)=\lambda G_e(t)E\left[S\right],$$(4) . For a non-empty system observed at time τ, we note a result by Mandjes and Uraniewski (2011) who analysed the approach to steady state identifying a function that behaves as the difference between the transient and stationary distributions in the limit, and of relevance in the following Section 3.2.

Example 3.2.

We use the Pareto distribution as it exhibits the heavy tailed non-exponential survival times observed in the prison domain, and we draw upon this later. For $G\sim \mathit{\text{Pa}}(\theta ,\alpha )$ with shift parameter $\theta >0$ and shape parameter $\alpha >0,$ for $x\ge 0,$ $G^c(x)={\theta }^{\alpha }{(x+\theta )}^{-\alpha }.$ The high variability of G is indicated by the fact that the tail decays as a power instead of exponentially. For $\alpha >1,\hspace{0.17em}E[S]=\theta {(\alpha -1)}^{-1},$ otherwise if $\alpha \le 1$ then S has infinite mean and $G^c(x)$ is not a directly integrable function. For $\alpha >2,\hspace{0.17em}\mathit{\text{Var}}[S]={\theta }^2\alpha {(\alpha -1)}^{-2}{(\alpha -2)}^{-1},$ for $1<\alpha \le 2$ the variance is infinite, and for α otherwise the variance is undefined. For $\alpha >1,$ $G_e^c(x)={\theta }^{\alpha -1}{(x+\theta )}^{1-\alpha }.$ If $\alpha >2$ then the SCV is $c_s^2=\alpha {(\alpha -2)}^{-1},$ then $E[S_e]=\theta {(\alpha -2)}^{-1}.$

3.2. Conditional distribution

Denote by τ the observation time, by $\tau +\delta$ the prediction time, and define the vector of the past arrival intensity, $\widehat{\lambda }=\{\lambda (t):0\le t\le \tau \}.$ For a stochastic process f(t) cut at an observation time τ, define the past process $\widehat{f}(t)=f(t)$ if $t<\tau$ and $\widehat{f}(t)=0$ if $t\ge \tau ,$ and define the future process as $\check{f}(t)=0,$ if $t<\tau$ and $\check{f}(t)=f(t)$ if $t\ge \tau .$ The past process can be thought of as an $M_t/G/\infty$ process in which arrivals are terminated at time τ, so we are able to draw on results by Goldberg and Whitt (2008) who, motivated by a model of a two-stage item inspection process, studied the distribution of the last departure time from a queue of a terminating arrival process.

For the analysis of the observed queue, we separate the contributions of the sub-populations of those observed at time τ and those arriving later. It is useful to consider the regions describing arrival and service pairs as depicted in Figure 2, where for each individual arriving to the system, a point is placed at (u, v), with the u-axis recording arrival times and v-axis recording service times. For example, Region (1)(1) $$G_e(t)=\frac{1}{E\left[S\right]}{\int }_0^t{G}^c(u)\mathit{\text{du}},$$(1) corresponds to individuals arriving and departing by τ, and the region $A_{(2\cup 3)}=\{(u,v)|\hspace{0.17em}u\le \tau ,u+v\ge \tau \}$ corresponds to individuals arriving on or before time τ. Denote by N(A) the number of arrival-service pairs in a region A. As per Theorem 3.1, the Poisson arrivals and general service times generate a Poisson process on the plane with the intensity of a point occurring at (u, v) being $(\lambda (u),g(v)).$ As independent splitting of Poisson processes results in Poisson processes, the numbers of points in two disjoint regions are independent Poisson random variables. For example, consider $A_{(3)}=\{(u,v)|\hspace{0.17em}u\le \tau ,u+v\ge \tau +\delta \},$ corresponding to individuals arriving by τ and present at $\tau +\delta .$ Then the proportion of individuals present at time τ whose remaining service time from that point is at least δ, is clearly $N(A_{(3)})/N(A_{(2\cup 3)}).$ This geometric depiction is captured in the following result.

Theorem 3.3.

(Goldberg & Whitt, 2008, Theorem 2.1) Conditional on there being n individuals in the system at time τ, the remaining service times are i.i.d., each distributed as a random variable $X_{\tau }$ with ccdf(5) $$G_{\tau }^c(x)=P\left(X_{\tau }>x\right)=\frac{1}{{\nu }_{\tau }}{\int }_{-\infty }^{\tau }\lambda (u)G^c(\tau +x-u)\mathit{\text{du}}.$$(5) where mean ${\nu }_{\tau }$ is given by(6) $${\nu }_{\tau }={\int }_{-\infty }^{\tau }\lambda (u)G^c(\tau -u)\mathit{\text{du}}={\int }_0^{\infty }\lambda (\tau -u)G^c(u)\mathit{\text{du}}.$$(6)

It directly follows that if $\lambda (t)=\lambda$ for $t\ge 0,$ and $\lambda (t)=0$ for t < 0, then ${\nu }_{\tau }=\lambda E[S]G_e(\tau )$ and (7) $$G_{\tau }^c(x)=\frac{G_e(\tau +x)-G_e(x)}{G_e(\tau )}.$$(7)

The above result provides an expression for the remaining service time distribution in Equation (5)(5) $$G_{\tau }^c(x)=P\left(X_{\tau }>x\right)=\frac{1}{{\nu }_{\tau }}{\int }_{-\infty }^{\tau }\lambda (u)G^c(\tau +x-u)\mathit{\text{du}}.$$(5) , which is required to calculate the conditional distribution in the following result.

Theorem 3.4.

(Weber, 2005, Theorem 6) The random variable $Q(\tau +\delta |\tau )$ with $Q(\tau )=n$ can be expressed as(8) $$Q(\tau +\delta |\tau )=\mathit{\text{Bi}}\left(n,p_{\tau }(\delta )\right)+\mathit{\text{Po}}\left(\check{m}(\tau +\delta )\right),$$(8) where(9) $$\check{m}(\tau +\delta )=E[\check{\lambda }(\tau +\delta -S_e)]E\left[S\right]={\int }_0^{\delta }\lambda (\tau +\delta -u)G^c(u)\mathit{\text{du}},$$(9) where $p_{\tau }(\delta )=G_{\tau }^c(\delta )$ is given in Equation (5)(5) $$G_{\tau }^c(x)=P\left(X_{\tau }>x\right)=\frac{1}{{\nu }_{\tau }}{\int }_{-\infty }^{\tau }\lambda (u)G^c(\tau +x-u)\mathit{\text{du}}.$$(5) , and $\mathit{\text{Po}}(·)$ and $\mathit{\text{Bi}}(·,·)$ are the Poisson and Binomial random variables, respectively.

The process $Q(\tau +\delta )$ can be written as $Q(\tau +\delta )=\widehat{Q}(\tau +\delta )+\check{Q}(\tau +\delta ),$ where $\widehat{Q}(t)$ has arrival rate $\widehat{\lambda }(t)$ and $\check{Q}(t)$ has arrival rate $\check{\lambda }(t).$ As constructed above, $Q(\tau +\delta |\tau )$ will be the number of points in regions $(3\cup 5)$ of Figure 2. $N(A_{(5)})$ has Poisson distribution with mean given in Theorem 3.1 and the result for $\check{Q}(\tau +\delta |\tau )$ is given by Equation (9)(9) $$\check{m}(\tau +\delta )=E[\check{\lambda }(\tau +\delta -S_e)]E\left[S\right]={\int }_0^{\delta }\lambda (\tau +\delta -u)G^c(u)\mathit{\text{du}},$$(9) . The distribution of $N(A_{(3)})$ is binomial with n trials and parameter $p_{\tau }(\delta )$ given in Equation (5)(5) $$G_{\tau }^c(x)=P\left(X_{\tau }>x\right)=\frac{1}{{\nu }_{\tau }}{\int }_{-\infty }^{\tau }\lambda (u)G^c(\tau +x-u)\mathit{\text{du}}.$$(5) .

Depending on the size of δ relative to τ, either existing or new arrivals can dominate the prediction. The calculations can reveal the contribution of each component to future demand requirements. For some quantities of interest the two sub-processes $Q(t)=(\widehat{Q}(t),\check{Q}(t)),$ from the expression in Equation (8)(8) $$Q(\tau +\delta |\tau )=\mathit{\text{Bi}}\left(n,p_{\tau }(\delta )\right)+\mathit{\text{Po}}\left(\check{m}(\tau +\delta )\right),$$(8) can be studied separately. The results extend easily to multiple classes, i, $\{Q_i(\tau )=n_i\}$ as each class is treated independently.

Proposition 3.5.

The conditional distribution is(10) $$P\left[Q(\tau +\delta )=y|Q(\tau )=n\right]=\sum _{k=0}^{\mathrm{\text{min}}\{n,y\}}P\left[\mathit{\text{Bi}}\left(n,p_{\tau }(\delta )\right)=k\left]P\right[\mathit{\text{Po}}\left(\check{m}(\tau +\delta )\right)=y-k\right],$$(10) with mean $E\left[Q(\tau +\delta )|Q(\tau )=n\right]=n{p}_{\tau }(\delta )+\check{m}(\tau +\delta )$, and variance $\mathit{\text{Var}}\left[Q(\tau +\delta )|Q(\tau )=n)\right]=n{p}_{\tau }(\delta )(1-p_{\tau }(\delta ))+\check{m}(\tau +\delta ).$

The departure process of the observed $M_t/G/\infty$ queue at time $\tau +\delta$ will be the superposition of a Poisson process with intensity ${\lambda }^{-}(\tau +\delta )={\int }_0^{\delta }\lambda (\tau +\delta -u)g(u)\mathit{\text{du}}$ and a binomial process with parameters n and $1-p_{\tau }(\delta )$ (Weber, 2005). A consequence of which is that the Poisson-in-Poisson-out feature of the unobserved $M_t/G/\infty$ system does not hold.

Under a high load, the binomial process will be the superposition of several point processes and will be approximately Poisson: for large n, (11) $$P[\tilde{Q}(\tau +\delta )=y|\tilde{Q}(\tau )=n]=P[\mathit{\text{Po}}(n{p}_{\tau }(\delta )+\check{m}(\tau +\delta ))=y].$$(11)

Thus considering the system began with an initial population drawn from a Poisson distribution with mean n and we employ this approximation in Section 4.

3.2.1. Last departure time

The behaviour of the decay of the initial population $\widehat{Q}(t)$ can be described by results of Goldberg and Whitt (2008). Suppose n individuals are observed at time τ and consider the time of the last departure. Define $M_n=\mathrm{\text{max}}\{X_{\tau }^1,\dots ,X_{\tau }^n\},$ where each $X_{\tau }^i$ is distributed according to Equation (5)(5) $$G_{\tau }^c(x)=P\left(X_{\tau }>x\right)=\frac{1}{{\nu }_{\tau }}{\int }_{-\infty }^{\tau }\lambda (u)G^c(\tau +x-u)\mathit{\text{du}}.$$(5) , then $P\left[M_n\le x\right]=G_{\tau }{(x)}^n.$ Now consider the case where the arrival process is terminated at time γ, but where an observation is not made. Following Goldberg and Whitt (2008), let D be the last departure time and define $T={(D-\gamma )}^{+},$ the remaining time after γ until the last departure. To determine the distribution of T, let $N=N_{\gamma }$ be a random variable with a Poisson distribution having mean ${\nu }_{\gamma },$ then $T\stackrel{d}{=}{M}_N.$

For any random variable Y with continuous cdf, let its quantile function be $q_Y\equiv q_Y(x)$ such that the number $P[T\le q_Y(x)]=x.$ Write $f(x)\sim g(x)$ as $x\to \infty$ when $f(x)/g(x)\to 1$ as $x\to \infty .$

Theorem 3.6.

(Goldberg & Whitt, 2008, Theorem 2.2) For any x > 0,(12) $$P\left[T\le x\right]=e^{-{\nu }_t{G}_{\gamma }^c(x)}.$$(12)

Then as $x\to \infty ,\hspace{0.17em}P[T>x]\sim {\nu }_{\gamma }{G}_{\gamma }^c(x),$ and for $e^{-{\nu }_{\gamma }}<x<1,$ (13) $$q_T(x)=q_{X_{\gamma }}\left(1-\frac{\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}(1/x)}{{\nu }_{\gamma }}\right).$$(13)

The moments of which can be calculated by $E[T^k]={\int }_0^{\infty }k{x}^{k-1}E[T>x]\mathit{\text{dx}}.$

Example 3.7.

For $G\sim \mathit{\text{Pa}}(\theta ,\alpha$), and with $\lambda (t)=\lambda$ for $t\in (-\infty ,\tau ],$ (14) $$q_T(x)=q_{X_{\tau }}\left(1+\frac{(\alpha -1)\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}(x)}{\lambda \theta }\right),$$(14) with (15) $$q_{X_{\tau }}(x)={\left(\frac{{\theta }^{\alpha -1}}{1-x}\right)}^{\frac{1}{\alpha -1}}-\theta .$$(15)

In Figure 9 for fixed $E[S]=3,$ we illustrate varying values $(\lambda ,c_s^2).$ The effect of increasing the SCV is seen in the quantiles of the last departure time T. These results can be employed in the application domain to study the decreasing content of a prison population with no new arrivals.

4. Empirical study and simulations

The statistical inference of infinite server queues has been studied given various types of incomplete data: (i) queue length data (Goldenshluger & Koops, 2019; Picklands, 1997); (ii) unmatched arrival and service times (Blanghaps et al., 2013); (iii) busy period data (Hall & Park, 2004); and (iv) arrival and departure counts (Li et al., 2019).

In this section, we present a method to generate predictions that consider the parameter uncertainty using Bayesian inference, considering queue length data from March 2015 to March 2019 on the Adult/Male sentenced population (MoJ, 2020b, Table 1.2b: Sentenced to immediate custody). We refer to the MoJ documentation for details on counting processes, but we interpret this as the sentenced population and only consider data from 2015 due to changes in the reporting processes.

We aim to provide a demonstration of the estimation approach for the model parameters using simulated monthly counts generated from published quarterly data, as the data is not sufficient to perform Bayesian inference about the system and model parameters. To generate monthly counts we: (i) divide the quarterly counts to obtain a monthly mean count; (ii) fit a smoothing spline to produce monthly predictions; (iii) add noise equal to the residual standard error from a linear model fit.

4.1. Bayesian analysis of the sentenced population

To specify the model, we adopt a linear arrival rate of the form $\lambda (t)={\beta }_0+{\beta }_{1}t$ and assume that the service time distribution for each offence group is Pareto, and can be adequately estimated from historical data (over a much longer time interval than the prediction lead time). Using public data (MoJ, 2020a, July 2018: June 2019) for the mean sentence length for each offence group, we specify that the service time is 50% of the sentence length for this demonstration.

For a large system, define the mean function $M(\tau +\delta )=E[\tilde{Q}(\tau +\delta )]$ using Equation (11)(11) $$P[\tilde{Q}(\tau +\delta )=y|\tilde{Q}(\tau )=n]=P[\mathit{\text{Po}}(n{p}_{\tau }(\delta )+\check{m}(\tau +\delta ))=y].$$(11) , and Pareto $(\theta ,\alpha )$ distribution as in Example 3.2, we have (16) $$v_{\tau }=\frac{({\beta }_0+{\beta }_1\tau )\theta }{\alpha -1}-\frac{{\beta }_1{\theta }^2}{(\alpha -1)(\alpha -2)},$$(16) (17) $$p_{\tau }(\delta )={\theta }^{\alpha -1}{(\delta +\theta )}^{1-\alpha }[1+\frac{{\beta }_1\delta }{({\beta }_0+{\beta }_1\tau )(2-\alpha )+{\beta }_1\theta }],$$(17) (18) $$\begin{array}{c}\check{m}(\tau +\delta )=\frac{({\beta }_0+{\beta }_1\tau )}{1-\alpha }[{\theta }^{\alpha }{(\delta +\theta )}^{1-\alpha }-\theta ]\\ +\frac{{\beta }_1}{(1-\alpha )(2-\alpha )}[{\theta }^{\alpha }{(\delta +\theta )}^{2-\alpha }-{\theta }^2-\delta \theta (2-\alpha )].\end{array}$$(18)

We specify prior distributions for arrival rate function and service distribution time parameters $\{{\beta }_0,{\beta }_1,\alpha ,\theta \},$ where arrival data is used to specify informative priors for arrival rate function. We define a Bayesian model for the monthly numbers of arrivals with non-informative priors. Denote by $Q_a(\tau +\delta )-Q_a(\tau )$ the number of individuals that arrive in $[\tau ,\tau +\delta$) and is Poisson distributed with mean ${\int }_{\tau }^{\tau +\delta }({\beta }_0+{\beta }_{1}u)\mathit{\text{du}}.$ We extract the posterior samples of β₀ and β₁. Since both posterior densities are bell-shaped, we define normal priors for β₀ and β₁ centred on the mean of posterior samples and standard deviation equal to the scaled posterior sample standard deviation for the model. We use a scaled posterior sample standard deviation in order to fully explore the parameter space during Markov chain Monte Carlo (MCMC). We specify a uniform prior for $\alpha \sim \text{Uniform}[2.5,10],$ where the lower bound avoids infinite variance and the upper bound ensures good convergence and mixing of the Markov chains. In practice, we observe that increasing the upper bound of α leads to non-identifiability. We use RStan to perform full Bayesian statistical inference by adopting the Hamiltonian Monte Carlo (HMC) and no-U-turn samplers (NUTS) (Carpenter et al., 2017).

We demonstrate long-term (k-step ahead prediction) and short-term (one-step ahead prediction) forecasts. The long-term forecasts are useful in informing long-term planning decisions in relation to the prison system capacity, whereas short-term forecasts could provide the insight needed for the day-to-day prison system operation. To generate predictions about the system behaviour, we perform the following steps:

1. Sample $({\beta }_0^{(k)},{\beta }_1^{(k)},{\alpha }^{(k)},{\theta }^{(k)})$ from the posterior distribution $[{\beta }_0,{\beta }_1,\alpha ,\theta |{\mathit{n}}^s]$ using RStan, where ${\mathit{n}}^s=(n_1,\dots ,n_{\tau })$ with $n_i,i=1,\dots ,\tau$ being the number of individuals present in the system at time i.

2. For each $k=1,\dots ,N$ compute $M^{(k)}(\tau +\delta )$ and sample ${\tilde{Q}}^{(k)}(\tau +\delta )$ from $\mathit{\text{Po}}(M^{(k)}(\tau +\delta )).$

3. Compute sample mean ${\mu }_{\tilde{Q}}(\tau +\delta )=\frac{1}{N}{\sum }_{k=1}^N{\tilde{Q}}^{(k)}(\tau +\delta ),$ which is our prediction.

4. Compute sample standard deviation, ${\sigma }_{\tilde{Q}}(\tau +\delta )=\sqrt{\frac{1}{N}{\sum }_{k=1}^N({\tilde{Q}}^{(k)}(\tau +\delta )-{\mu }_{\tilde{Q}}(\tau +\delta ){)}^2},$ which represents uncertainty about our prediction.

4.2. Predictions and simulation results

Example 4.1.

We demonstrate the Bayesian model specification to predict the number of prisoners in the Theft offence group where $E[S]=5.22$ months. From the arrival count data, we obtain the posterior sample for β₀ with mean 1376.5 and standard deviation 9.7. Similarly, the mean and standard deviation of posterior sample for β₁ are -11.5 and 0.3 respectively. We then specify the priors: ${\beta }_0\sim \text{Normal}(1376.5,10\times 9.7),\hspace{0.17em}{\beta }_1\sim \text{Normal}(-11.5,10\times 0.3)$ and $\alpha \sim \text{Uniform}[2.5,10],$ with $\theta =5.22(\alpha -1).$ We set two Markov chains with 10,000 iterations for each chain (including warmup).

In Figure 3, we produce the density plots of marginal posterior distributions for model parameters. We observe that the posterior sample density of β₀ is bell shaped and centered around 677.4 and the standard deviation of posterior samples is 17.76. Similarly, the values of β₁ is centred around −3.77 with posterior standard deviation 0.54. We observe that under our model specification the number of arrivals in Theft Offence group gradually decreases over time. The distributions of posterior samples for θ and α are less informative, since we included the mean of service time in our prior specification.

In Figure 4 we illustrate long and short term predictions from August 2019 to March 2020. We note that the short-term projections are more computationally expensive as in order to obtain a new observation, we update the posterior distribution by re-running a Stan programme. The left panel plot in Figure 4 shows the long-term forecast (for 8 months). The number of offenders in Theft offence category is expected to decrease over time. Our predictions slightly underestimate the true values, however the true values lie within two standard deviation prediction interval. We observe that uncertainty about our projections increases with time. The right panel plot in Figure 4 shows the short-term forecasts. The predictions are closer to the true values together with the slightly lower uncertainty about our predictions. To access the performance of the proposed model, we use the Root Mean Square Error (RMSE): $\mathit{\text{RMSE}}=\sqrt{q^{-1}{\sum }_{\delta =1}^q(n_{\tau +\delta }-{\mu }_{\tilde{Q}}(\tau +\delta ){)}^2},$ where a lower value indicates a better model performance. The RMSE for long-term predictions is 53.08, whereas for short-term predictions is 20.9.

Example 4.2.

Adopting a similar approach as in Example 4.1, we consider the number of offenders for two further offence groups: Sexual offences and Public order offences. The predictions for Public order offences are illustrated in Figure 5. For Sexual offences the RMSE for long-term forecast is 7.16, whereas for short-term forecast is 5.10. For Public order offences: the RMSE for long-term forecast is 2.97, whereas for short-term forecast is 1.89.

In the next example we reuse the posterior samples of model parameters derived in Examples 4.1 and 4.2 to demonstrate how the model can provide insight into policy modifications.

Example 4.3.

A change point can be studied, representing switching from one kind of service to another at time τ (that does not involve an observation), the model is an $M_t/G^o,G^{\nu }/\infty ,$ (Korhan Aras et al., 2017) (using the superscript o for old and ν for new). The distribution of $Q(\tau +\delta )$ is Poisson with mean $\widehat{m}(\tau )G_e^{o,c}+\check{m}(\tau +\delta ).$

In Figures 6, 7 and 8 we simulate the effect of increasing and decreasing the mean service time on the offence group population for theft offences, sexual offences and public order offences, respectively, and where the solid lines are the predictors and the dashed lines represent two standard deviation prediction intervals. We remark that to produce these results we reuse samples from the posterior distributions of ${\beta }_0,{\beta }_1$ and α. The new values of θ are given by $\theta =E[S^{\nu }](\alpha -1).$

4.3. Discussion

In the previous section we illustrated how to use the Bayesian model for short and longer-term predictions, and to provide insight into the implications of possible policy modifications. In Appendix B, drawing on other theoretical results, we briefly describe how modification of sentencing and custody rates can be seen as a form of intervention to enable congestion event recovery. Empirical demonstration of the value of such analysis is beyond the scope of the paper, as this would require data not available to us.

With regard to the predictions in Section 4.2, we note that time-series forecasting methods such as ARIMA models (Shumway et al., 2000) can predict the future prison population by describing the autocorrelation in the data. A direct comparison, using the data from Section 4.1, is presented in Appendix C, and in Section 5 we provide some comments on the different approaches to forecasting.

We now briefly discuss the impact of the assumptions of the mathematical model described in Section 3. It is assumed that the time served by individuals at the observation point t = 0 is unknown, but if the elapsed service times $\{y_i:i=1,\dots ,n\}$ of the observed population $\widehat{Q}(t)$ are recorded, it can be more effective to use this information to make predictions, while carrying a greater cost (Duffield & Whitt, 1997). In contrast to Equation (5)(5) $$G_{\tau }^c(x)=P\left(X_{\tau }>x\right)=\frac{1}{{\nu }_{\tau }}{\int }_{-\infty }^{\tau }\lambda (u)G^c(\tau +x-u)\mathit{\text{du}}.$$(5) , the conditional remaining service time ccdf for elapsed service time x is $H_x^c(t)=G^c(t+x)/G^c(x).$

As the conditional remaining service times are no longer identically distributed, the complication of the distribution $\widehat{Q}(t)$ increases and the resulting process is a Markov process (Korhan Aras et al., 2017), although the estimates for the mean number remaining at time t and variance have simple forms: $E[\widehat{Q}(t)]={\sum }_{i=1}^n{H}_{y_i}^c(t),\mathit{\text{Var}}[\widehat{Q}(t)]={\sum }_{i=1}^n{H}_{y_i}^c(t)H_{y_i}(t).$

Further, as noted by Whitt (1999), the importance of conditioning upon the time served at an observation point increases as the service-time distribution differs more from an exponential distribution. If G is highly variable, then the elapsed holding time can greatly help in future prediction and a long elapsed holding time tends to imply a very long remaining holding time. Let $Y(\alpha ,\theta )$ denote the Pareto service time distribution G as defined in Example 3.2 and let $Y_x(\alpha ,\theta )$ the associated ccdf for elapsed service time x. By the result for $H_x^c(t)$ above, Duffield and Whitt (1997, Theorem 8) showed that $H_x^c(t)=G^c(t{(\frac{x}{\theta }+1)}^{-1}),$ which implied $Y_x(\theta ,\alpha )\stackrel{d}{=}\left(1+\frac{x}{\theta }\right)Y(\theta ,\alpha ).$ Hence $E[Y_x(\theta ,\alpha )]=\left(1+\frac{x}{\theta }\right)E[Y(\theta ,\alpha )],$ that is, the mean remaining service time $E[Y_t(\alpha ,\theta )]$ is approximately proportional to the elapsed service time t.

5. Concluding remarks

This work was motivated by the problem of predicting the population of housed inmates within the prison system in England and Wales, where the size of the prison population is recorded on a regular basis, with attention both to short-term predictions to enable local resource planning, and longer-term projections that can provide insight to policy makers regarding the impact of potential policy variations.

We studied the transient behaviour of the time-varying infinite server queue, $M_t/G/\infty ,$ fed by a non-homogeneous Poisson arrival process whose occupancy is observed at discrete points in time, but the time in service to that point is not known. Drawing on this analysis, and using publicly available data, we built a model that could be used as a decision support tool for custodial elements of the prison system. We illustrated the use of such a model for population prediction and for analysing the implications of changing external factors influencing the prison population such as government guidelines and sentencing policies. The proposed model produces predictions together with uncertainty bands and aligns with the current guidelines on informed decision making in the UK government (Treasury, 2015).

It is beyond our scope to compare the queueing theory approach to generating predictions with multiple forecasting methods, but we note recent work arguing that for time-series methods Liu et al. (2021) to support general scenario analysis requires extracting components known as “features” from the time-series data, that in turn can be used to generate alternative scenarios (Kegel et al., 2017; Tuominen et al., 2022). In the queueing theory approach, model parameters have a direct interpretation for the application domain, which straightforwardly enables the study of a variety of public policy initiatives by setting different input values to the model. We demonstrated this in Example 4.3 for changes in the distribution of sentence lengths. As a further example, changes under consideration to the sentencing and release of serious and dangerous sexual and violent offenders (MoJ, 2019), could be studied by changes in the arrival rate (e.g., the custody rate, conviction rate), studies we were unable to undertake, as they require data not available to us. In contrast, time-series methods are not so amenable to what-if style scenario analysis as they offer the possibility of forecasting future observations, but with limited interpretability of the fitted model (Petris et al., 2009).

The contributions of this paper are: (i) the novel synthesis of results from several authors about transient and stationary behaviour of the $M_t/G/\infty$ queue to enable construction of a model suited to short and longer-term predictions, and to supporting considerations of parameter uncertainty; and (ii) illustration of the approach to potential policy changes in the real-world domain of prison occupancy.

Reflecting the data available for model building, we focused on the situation where the system has non-empty initial state and where the elapsed time of each individual in the system is not known. The dynamics of the queue is a combination of those already in service at some time and those who subsequently arrive, and separation into initial content and new input allowed these sub-populations to be analysed jointly and separately. We drew on results for the transient and stationary distributions of these queueing systems from several authors to enable an analytic approach. Then, using a Bayesian approach with public historical data that allows the inclusion of expert knowledge, we considered parameter uncertainty involved in the prediction of future arrivals and presented a model that maintains interpretability for the domain application. Incorporating other sources of uncertainty into our process Whitt (2002), including model and process uncertainty, and quantifying the contribution from each within our application is a topic of future research.

Further, we note that restoration of the departure process as approximately Poisson also allows the approach to be extended to a network of processes, referred to as a ${(M_t/G/\infty )}^N/M$ model, in which queue length distribution models have time-dependent product form and would be appropriate for models of many service systems (Massey & Whitt, 1993), including further development of a model specific to the prison domain.

The approach is potentially applicable to other service systems, but the queueing model properties of interest will vary according to the application context, and the choice of what approach to take will also depend on the available data. Data driven development of a system model is an increasingly popular approach as discussed in Mandelbaum et al. (2019). However, the infrastructure to collect and manage data at multiple phases and timescales is not yet comprehensive, so parameter inference remains a challenging problem because of limited data availability about successive system states, and for model building methods that can take advantage of available but incomplete data are essential. We acknowledge that some of our articulated modelling assumptions may not apply in domains where fine-grained staffing implications are of interest, for example, hospitals and call centres.

The availability of future information (via predictive algorithms, machine learning, or local observation) and the resultant novel queueing analysis, has policy design implications, described as “a broader shift from being reactive to proactive” as future information becomes part of the policy maker’s toolkit (Spencer et al., 2014; Walton & Xu, 2021).

Acknowledgements

We would like to thank our colleagues from the Ministry of Justice for helpful discussions in the development of this work, and the journal reviewers for constructive comments that helped to improve the presentation. The work was conducted as part of the Managing Uncertainty in Government Modelling project supported by The Alan Turing Institute. This work was supported by the Additional Funding Programme for Mathematical Sciences, delivered by EPSRC (EP/V521917/1) and the Heilbronn Institute for Mathematical Research.

Disclosure statement

No potential conflict of interest was reported by the authors.

Appendix A:

Prison occupancy

Phases of the prison system

Following an arrest by the police and being charged with a crime, individuals are classified within 12 offence main groups (Table C1), which naturally form different service classes. The police decide if a person can be released on bail, but may have to follow certain rules until a hearing at a court, or if they are held on remand in prison until a hearing at a court. If an individual is found guilty by a court, a custodial sentence can be imposed. A sentence can be suspended, or served in the community and the vast majority of sentenced offenders do not spend their whole sentence in custody. The sentenced population is made up of Determinate sentenced and Indeterminate sentenced. Each sentence is marked as one of two types: a standard determinate sentence (SDS) or an extended determinate sentence (EDS). Most offenders are given an SDS, which means that half way through their sentence they are released from prison and serve the remaining time in the community. An EDS means an offender serves 2/3 of their sentence in prison before being considered for release. Sentenced offenders who have been released from prison before the end of their sentence can be recalled to custody for breach of release conditions.

Figure C1 displays a simplified version of prisoner journeys through the custodial system, with the total prison population divided into three main holding phases (displayed percentages are as at June 2019) (MoJ, 2019): (i) on remand (11%), (ii) sentenced prisoners (79%) and (iii) on recall (9%). Prisoners within the licence phase are in the community.

Arrivals to the Remand phase, are recorded in two streams, Untried prisoners and Sentenced prisoners. Individuals remain in remand for some time according to a length of stay distribution $G_1.$ Upon completion of the Remand phase, the model assumes all Sentenced prisoners proceed to the Sentenced population and a proportion of the Untried prisoners, with the remaining released. An arrival stream from the courts ${\lambda }_2(t)\ge 0$ brings individuals directly into the Sentenced population. All departures from the Sentenced population are released into the community, and are all under Licence, for some time, during which they can be Recalled into custody. If an individual is recalled, they re-enter the custodial system and this figure assumes a prisoner can only be recalled once. In practice, this is not the case, although the data available does not provide information on repeatedly recalled prisoners.

Data

Individual prison establishments record details of individuals on a prison IT system and this data is collected as a central database managed by MoJ and the HM Prison and Probation Service (HMPPS) to produce the various analyses of prison population. Statistics are regularly released as well as projections of the prison population (Ministry of Justice, 2016, 2018).

Between 1993 and 2012 the prison population in England and Wales increased by 41,800 prisoners to over 86,000, studies of this growth can point to patterns in the data that align with earlier identifiable policy changes. Such examples illustrate how modelling approaches enabling an understanding of the impact of possible policy changes can be valuable in this domain.

Since 1993 the courts sentenced more offenders each year. Almost all of the increase (98%) took place within two segments of the population - those sentenced to immediate custody (85% of the increase) and those recalled to prison for breaking the conditions of their release (13% of the increase). Three offence groups, violence against the person, drug offences and sexual offences have had a particular impact on the prison population: the number of these offenders increased and the mean length of stay increased.

The proportion of offenders whose sentences were 4 years or more, grew to 54% of the sentenced prison population in 2012, compared to 45% in 1993. The mean sentence length over all offences in 2018 was 17.3 months. For more serious offences, the mean prison sentence was 58.3 months, more than 24 months longer than in 2006. More than two and a half times as many people were sentenced to 10 years or more in 2018 than in 2006.

From 1999 to 2011, the mean time served increased from 8.1 to 9.5 months for those released from determinate sentences. This was likely due to an increase in the mean custodial determinate sentence length handed down by the courts, and a decline in the parole release rate which meant that offenders had served longer by the time they were released. The second largest increase was within the recall population, caused by changes to the law making it easier to recall prisoners, and changes introduced in the Criminal Justice Act 2003 which lengthened the licence period for most offenders.

Appendix B:

Steady state recovery from congestion

To illustrate the applicability of further theory to support policy analysis, we present Duffield and Whitt (1997)’s analysis of what is termed congestion events, that is, where an unusually high number of individuals has been observed, and allowing one to investigate the effects of possible interventions (such as change of arrival rate, or pausing arrivals) on recovery from congestion. These results rely on an assumption that the system is in steady state, i.e., $\lambda (t)=\lambda ,$ which could be motivated by a slowly changing arrival rate, or for scenario analysis. These results follow directly from the results for the model described in Section 3.2 and can be employed to characterise the effect of possible interventions that could be introduced to change the rate of recovery.

For all $\tau ,$ $\nu ={\nu }_{\tau }=\lambda E[S],$ $\check{m}(\tau +\delta )=\lambda E[S]G_e(\delta )$ and $p_{\tau }(\delta )=G_{\tau }^c(\delta )=G_e^c(\delta )$ and the conditional mean is (B1) $$E\left[Q(\tau +\delta )|Q(\tau )=n\right]=\lambda E\left[S\right]G_e(\delta )+n{G}_e^c(\delta ).$$(B1)

As Goldberg and Whitt (2008) noted, since ${\nu }_{\tau }$ and $G_{\tau }^c(\delta )$ can be computed, it can be directly determined when the steady state approximation is reasonable.

Duffield and Whitt (1997) proposed the steady state rare event, as an approximation for what is seen at a random (steady state) time and also as an approximation for the associated hitting time event. The likelihood of a high congestion event, that is, having a large number N of individuals present, $P\left[N\ge n\right]={\sum }_{k=n}^{\infty }P[\mathit{\text{Po}}(\nu )=k]$ with $\nu =\lambda E[S]$ at $t=0.$

Define a recovery level k, then a system has recovered from a rare congestion event at time 0 when $Q(t)$ first reaches k with $\nu <k<n.$ To define the recovery time for the mean, let (B2) $$\beta \equiv {\beta }_{n,k}=\mathrm{\text{inf}}\left\{t\ge 0:E\left[Q(t)|Q(0)=n\right]\le k\right\}.$$(B2)

The following result describes how the system recovers from congestion and the way it depends on cdf G.

Theorem B.1.

(Duffield & Whitt, 1997, Theorem 2) If there is no t such that $G^c(t^{-})>0=G^c(t)$, then $G_e^c(t)$ is continuous and strictly increasing. Then the recovery time for the mean, $\beta$ in Equation (B2)(B2) $$\beta \equiv {\beta }_{n,k}=\mathrm{\text{inf}}\left\{t\ge 0:E\left[Q(t)|Q(0)=n\right]\le k\right\}.$$(B2) is the unique root of the recovery equation,(B3) $$G_e^c(t)=\frac{k-\nu }{n-\nu }.$$(B3)

More generally, $\beta =\mathrm{\text{sup}}\left\{t:G_e^c(t)>\frac{k-\nu }{n-\nu }\right\}.$

This result is obtained by rearranging Equation (B1)(B1) $$E\left[Q(\tau +\delta )|Q(\tau )=n\right]=\lambda E\left[S\right]G_e(\delta )+n{G}_e^c(\delta ).$$(B1) . It can be seen that the impact of G on $\nu$ depends on whether recovery level k is close to n or close to v. When the recovery level k is close to v, the cdf matters greatly. If the recovery level k, is closer to n than v, then the recovery time v depends more on the initial portion of the cdf $G_e$ than upon its tail. Consistent with this observation, the impact of a long-tail cdf on v differs little from the impact of a short-tail cdf with the same mean if k is suitably close to n, but it differs dramatically if k is suitably close to v. However, “suitably close” depends on the decay rate of the ccdf $G_e^c.$

If the system does experience a congestion event, the recovery can be aided by reducing the arrival rate, or by pausing arrivals until the congestion has cleared. If the arrival rate is reduced to a new constant rate, is equivalent to reducing the mean, as $\nu =\lambda E[S],$ produces a new recovery equation of the form in Equation (B3)(B3) $$G_e^c(t)=\frac{k-\nu }{n-\nu }.$$(B3) . For $k<m<n,$ clearly $(k-\nu )/(n-\nu )$ is increasing in $\nu ,$ so that the recovery time for the mean is reduced by decreasing $\nu .$ If the arrivals are terminated until recovery is achieved, $\nu =0$ in the recovery equation. Then further recovery can be achieved by turning the arrivals on and replacing $n$ with $k,$ and $j$ for $k$ in the recovery equation, such that $m<j<k.$

Example B.2.

For $G\sim \mathit{\text{Pa}}(\theta ,\alpha$), the recovery equation in Equation (B3)(B3) $$G_e^c(t)=\frac{k-\nu }{n-\nu }.$$(B3) we have $t={(\frac{{\theta }^{\alpha -1}(n-\nu )}{k-\nu })}^{\frac{1}{\alpha -1}}-\theta .$ In Figure C10, for fixed $E[S]=3$ and recovery level $k=\lceil \nu +1\rceil ,$ the black lines illustrate the recovery time for varying levels of congestion for different scenarios $(\theta ,\alpha ,c_s^2),$ where the x-axis shows the congestion level relative to the steady state mean $\nu .$ To aid the recovery for the three scenarios (black lines), we reduce the arrival rate by 20 percent, that is, ${\lambda }_2=0.8{\lambda }_1$ and illustrate the effect on the corresponding dashed lines (blue lines).

Appendix C:

Comparison to time series analysis

We apply traditional time-series methods, ARIMA models (Shumway et al., 2000) to produce prison population predictions for the offence groups considered in Section 4.2. The order and the values of parameters are obtained using R package forecast (Hyndman and Khandakar, 2008). The forecasts of the number of prisoners in the Theft offence group are produced by ARIMA(2,1,0) with drift given by: $$y_t=-9.76+0.94y_{t-1}+0.462y_{t-2}-0.41y_{t-3}+{ϵ}_t.$$

Similar to the analysis performed in Section 4.2, one-step ahead predictions are produced by iteratively adding one data point at a time and refitting the model and we note that the order of the model was updating with the changes in the training data set. The RMSE for long-term predictions is 30.21, whereas for short-term predictions is 10.16. The results are presented in Figure C11, and we observe that the length of the prediction intervals is lower compared to the predictions produced in Figure C4. This is due to the fact that the parameters of ARIMA model are set to the fixed quantities, whereas we assign probability distributions for the parameters in our proposed model. Similar results are presented for the Sexual offences and Public order offences groups in Figure C12. For example we present the $k$-step ahead forecasts for the Sexual offences group which are produced by ARIMA(2, 2, 2): $$y_t=2.55y_{t-1}-1.95y_{t-2}+0.23y_{t-3}+0.16y_{t-4}+{ϵ}_t-1.62{ϵ}_{t-1}+0.68{ϵ}_{t-2},$$ where the RMSE for long-term predictions is 4.71, whereas for short-term predictions is 4.42, where the above comment holds observing lower prediction intervals. The $k$-step ahead forecasts for the Public order offence group was produced by ARIMA(2,0,2): $$y_t=93.22+0.28y_{t-1}+0.47y_{t-2}+{ϵ}_t+0.91{ϵ}_{t-1}+0.85{ϵ}_{t-2},$$ where the RMSE for long-term predictions is 5.44, whereas for short-term predictions is 2.16.

Figure 1. Phases of the prison system.

Figure 2. Regions describing arrival and service pairs (u, v).

Figure 3. Left to right: density plots of marginal posterior distributions for parameters β₀, β₁, α and θ.

Figure 4. Prediction of the number of prisoners in the Theft offences group from August 2019 to March 2020. Left panel: long-term projections (red line) together with two standard deviation prediction interval (grey shaded region). Right panel: short-term projections (red line) together with two standard deviation prediction interval (grey shaded region). The observed values are in black.

Figure 5. Prediction of the number of prisoners in the Sexual offences (top row) and Public order offences (bottom row) groups from August 2019 to March 2020. Left panel: long-term projections (red line) together with two standard deviation prediction interval (grey shaded region). Right panel: short-term projections (red line) together with two standard deviation prediction interval (grey shaded region). The observed values are in black.

Figure 6. Example 4.3. Sentencing policy modification for Theft offences group where $E[S^o]=5.22.$ The red line is associated with unchanged parameters, the blue line is associated with $E[S^{\nu }]=8$ and the green line with $E[S^{\nu }]=3.$ Solid line is the predictor. Dashed lines represent two standard deviation prediction intervals.

Figure 7. Example 4.3. Sentencing policy modification for Sexual offences group where $E[S^o]=29.71.$ The red line is associated with unchanged parameters, the blue line is associated with $E[S^{\nu }]=48$ and the green line with $E[S^{\nu }]=24.$ Solid line is the predictor. Dashed lines represent two standard deviation prediction intervals.

Figure 8. Example 4.3. Sentencing policy modification for Public order offences group where $E[S^o]=3.79.$ The red line is associated with unchanged parameters, the blue line is associated with $E[S^{\nu }]=12$ and the green line with $E[S^{\nu }]=1.$ Solid line is the predictor. Dashed lines represent two standard deviation prediction intervals.

Figure 9. Example 3.7. For fixed $E[S]=3,$ varying $(\lambda ,c_s^2).$

Figure C11. ARIMA model prediction of the number of prisoners in the Theft offences group from August 2019 to March 2020. Left panel: long-term projections (red line) together with two standard deviation prediction interval (grey shaded region). Right panel: short-term projections (red line) together with two standard deviation prediction interval (grey shaded region). The observed values are in black.

Figure C10. Example B.2. Congestion level relative to the mean $n\nu$ vs. recovery time $t.$ For fixed $E[S]=3,$ $\nu =30$ and varying $(\theta ,\alpha ,c_s^2)$ (black lines). An intervention of ${\lambda }_2=0.8{\lambda }_1$ is illustrated by the corresponding blue lines.

Figure C12. ARIMA model prediction of the number of prisoners in the Sexual offences (top row) and Public order offences (bottom row) groups from August 2019 to March 2020. Left panel: long-term projections (red line) together with two standard deviation prediction interval (grey shaded region). Right panel: short-term projections (red line) together with two standard deviation prediction interval (grey shaded region). The observed values are in black.

Table 1. Offence groups.

1.	Violence against the person	7.	Possession of weapons
2.	Sexual offences	8.	Public order offences
3.	Robbery	9.	Miscellaneous crimes against society
4.	Theft Offences	10.	Fraud Offences
5.	Criminal damage and arson	11.	Summary Non-Motoring
6.	Drug offences	12.	Summary motoring