1 Introduction

Public transport system is the mode of travel for the common man. With the urbanization and population explosion, traffic management is becoming a challenging task. It undergoes a lot of changes in the recent times towards automation and intelligent fleet management. Artificial Intelligence and Machine learning based systems will revolutionize future transportation. One key aspect of automation is the prediction of passenger traffic to reduce congestions, accident rates, waiting time and to increase passenger comfort. Accurate passenger prediction helps in proper bus and crew scheduling. Thus the constrained services can be effectively managed, thereby improving the quality of the service.

The real time management of bus services enhances the profitability of the service providers and reduces the waiting time of passengers. The more efficient the public transport, the more the common man gets attracted towards it, reducing the usage of private cars, traffic congestion and pollution. The passenger traffic is a non stationary stochastic process. The morning traffic is mostly contributed by working persons and the student community. The noon traffic includes the common man who is travelling for their daily needs and the post noon session covers all those who are returning home after their work. The bus accessibility, type of traffic, land use, trip timing are all the influential factors affecting passenger flow. The analysis of public transport systems include geospatial data analysis, ticket data analysis and crew analysis. Often, such data are stochastic which requires statistical modelling and characterization.

In the present study, a set of Pune Mahanagar Parivahan Mahamandal Limited (PMPML) bus routes are analyzed with the help of Gaussian Process Regression (GPR) with the objective to maximize the earning per kilometer (epk) and passenger comfort. The GPR is a powerful class of Machine learning algorithm that relies on a few parameters to make predictions making it an ideal candidate for predicting the passenger traffic with limited data.

2 Literature Review

A set of models including Auto Regressive Moving Average models (ARMA), error correction models (VECM and EC-VARMA), space time ARMA (STARMA), and multivariate autoregressive space state models (MARSS) are used in the United States [1] for traffic analysis and prediction. GPS data collected from private cars is used to model a multivariate state space ARIMA model [2, 3] to forecast travel times in many countries. Short term passenger data prediction, [4] based on time series data, achieved better prediction performance. Neural networks [5] are employed to predict traffic data with the trend towards intelligent systems.

GPR has been found to be efficient on several multi-dimensional regression tasks [6]. GPR is one of the new machine learning method in the field of Bayesian Statistics [7]. Gaussian Processes (GPs) have been applied successfully to the problem of learning predictive state models [8]. GPs are now being used in modelling traffic data and prediction [9] as the traffic pattern exhibits a random behaviour. The wide range of applications in which GPR is used includes controlling MLS of Maglev train [10], rainfall prediction [11], solar PV panel power prediction [12], predicting covid-19 pandemic confirmed cases in South Korea [13] and assessing the energy behaviour of RC shear walls [14]. Several technologies are being incorporated to improve data collection [15], energy management [16] and traffic efficiency [17]. The Gaussian process regression technique is now widely used in many areas like rural traffic prediction [18], traffic modelling in high speed networks [19], Stock price forecasting [20], energy measurement and verification [7], soil moisture modelling [21], predicting multi-step ahead SPEI drought index [22], weather forecasting [23] and even in big data analysis [24].

The Gaussian process regression is currently being used in the prediction of both passenger and data traffic. The GPR technique is used in air traffic for quantifying data driven trajectory accuracy and uncertainty in Flight trajectory predictions [25, 26]. In maritime traffic, GPR is used for traffic probabilistic prediction based on ship motion pattern extraction [27] that includes both destination and trajectory prediction. In rail transport, the average speed is one of the most critical factors affecting travel time between two access points. GPR technique is used for the estimation of the average speed for a railway signalling system [28]. Considering the various aspects of traffic management, the global speed prediction [29] and individual longitudinal driving behaviour description of drivers [30], based on historical data, is done using multiple Gaussian process regression for different driving styles.

3 Theory

GPR is a supervised machine learning technique that provides a mapping from input to output[31]. The GPR techniques are more application specific than general regression models. The training procedure adapted by GPR process is that it considers the whole training data set at each time the model makes a prediction. Thus the adversial problems are completely eliminated when compared with traditional machine learning approaches. The response of GPR model is robust against slight variations in the passenger data. The main benefit of GPR is that it works well with small datasets and provide uncertainty measurements on the predictions. The training process using GPR model is relatively simpler compared to other ML models, since the hyperparameters of GPR are adaptively obtained by maximum likelihood estimation process. The GPR based models can easily overcome the overlearning or other data driven issues, since the training process considers the whole set of training data. The predictive accuracy, that is quantitatively constructed is least affected by any Gaussian noise [32]. Thus GPR is selected for the prediction of passenger traffic.

3.1 Passenger Arrival Model

It is required to maintain uniform number of passengers in all stages to increase the passenger comfort and revenue. The passenger arrival is modelled to accomplish this, so that proper bus scheduling is done. The arrival of passengers is conjectured as a sum [33] of several time varying Poisson process [3], where the probability of x passengers in a given period given by

$$\begin{aligned} P(x,\lambda _i)&=\sum _{i=1}^{n}\frac{ e^{-\lambda _i}\lambda _i^{x}}{x!} \end{aligned}$$
(1)

Here \(\lambda (t)\) is the average passengers per interval which is a random process. The passengers arriving at various stops choose buses for specific routes accordingly, that account for the bus passenger traffic, whose stochastic nature is analyzed below.

3.2 Passenger Traffic Model

The passenger traffic at any instant shows random behaviour and does not depend on the previous passenger history, making it a Markovian process. Linear estimation models like Kalman filter [34], Bayesian network based filters like particle filter [35] are used to model passenger traffic. The theory behind these bayesian filters (Kalman and Particle) is the estimation of posterior distribution of the passenger data by iteratively computing over the state of a dynamical system as

$$\begin{aligned} p(y_k|z_{1:k},u_{1:k-1})\propto p(z_k|y_k)\times \int p(y_k|x_{k_1},u_{k_1})p(y_{k_1}|z_{1:k})\,dy_{k1} \end{aligned}$$
(2)

Here y represents the true state and z represents the measurements of a Hidden Markov Model (HMM). The true state y is an unobserved Markov process obtained from the passenger data, and the measurements z are observations of HMM, represent the passenger count. The parameter \(u_{1:k-1}\) are the initial control measures or hyperparameters. The prediction model described by Eq. 2, \(p(z_k|y_k)\) represents the likelihood of the passenger data and \(\int p(y_k|y_{k1},u_{k1})p(y_{k1}|z_{1:k})\,dy_{k1}\) indicates the prior, whose product results in the posterior distribution. Generally, these filter models are parametric representation of the processes involved [36]. The major challenge that is faced in parametric models is that, they are not intelligent enough to reproduce all the aspects of a dynamical system. So the need of the hour is a non-parametric Bayesian model. Gaussian process regression(GPR) overthrows these parametric models and learns the prediction and observation models dynamically.

3.3 Gaussian Process Regression

Gaussian processes are non-parametric models that define a distribution over the mapping function f(x). The mapping function is a random variable of the input passenger data \(\textit{X}\). The distribution [37] is over infinite dimensions, but the computations are done using limited finite resources. This property makes GPs suitable for many spatial and temporal prediction tasks. Let x be the variable associated with the mapping function and \(x'\) be its time shifted version. A GP prior on the mapping function f(x) is completely specified by its mean and covariance as

$$\begin{aligned} f(x)\sim GP(m(x),k(x,x')) \end{aligned}$$
(3)

Here GP represents a Gaussian process distribution with mean m(x) and \(k(x,x')\) being the kernel or covariance function. \(k(x,x')\) gives the covariance of f(x) and \(f(x')\), and k(x, x) gives the variance of f(x). The regression task aims at predicting the target \(f(x_*)\) for new passenger data \(x_*\), that does not belong to the input data set \(\textit{D}\) (where \(\textit{D}\) \(\in\) \(\textit{X}\)). To predict this value, GPR model uses the GP prior on the mapping function f(x), along with the knowledge provided by the data set \(\textit{D}\) to calculate the posterior distribution over the mapping function. Finally, it makes inferences about the conditional distribution of the function value at \(x_*\).

The generalization properties of GPs rest almost entirely within the choice of kernels. The kernel specifications are computed using the passenger dataset and not inbuilt ones. The passenger data, modelled as a discernible time series, is used for modelling the prior which in turn is employed for computing the kernel parameters. The joint distribution of f(x) and \(f(x_*)\) is multivariate normal represented as

$$\begin{aligned} \left[ \begin{aligned} f(x) \\ f(x_*) \end{aligned}\right] \sim N \left( \left[ \begin{aligned} m(x)\\ m(x_*) \end{aligned}\right] ,\left[ \begin{aligned} k(x,x') k(x_*,x) \\ k(x_*,x) k(x_*,x_{*}^{'} \end{aligned}\right] \right) \end{aligned}$$
(4)

The conditional distribution is then calculated for predicting the function values \(f(x_*)\) at new test points \(x_*\), given by

$$\begin{aligned} \begin{aligned} f(x_*)|f(x)\sim N(k(x_*,x)k(x,x)^{-1}|f(x)&-m(x)|+m(x_*),k(x_*,x_*)\\ {}&-k(x,x_*)k(x,x)^{-1}k(x,x*)) \end{aligned} \end{aligned}$$
(5)

The methodology for analysing the bus passenger traffic using the above model is detailed in the next section.

4 Methodology

The methodology of work is as shown in Fig. 1. The different steps are

  • Data mining and modelling

  • Modelling of prior distribution

  • Gaussian process modelling

  • Comparison with existing methods

  • Interpretation and analysis of results

Figure 1
figure 1

Methodology of work.

Full size image

4.1 Data Mining and Modelling

The data related to the Pune city bus traffic is provided by the PMPML authorities in various formats. The data with their specifications from various stakeholders in the organization is given in Table 1.

Table 1 Data collected from PMPML.
Full size table

The Geospatial data gives the latitude and longitude values of different stops in the Pune city for the selected routes. The tripsheet report gives the information regarding the number of trips assigned in each route, direction of travel (up or down trip), the allotted bus number and crew members in each trip, trip timings and the assigned depot name. The passenger profile reports include the details of the number of passengers (including children) boarded from different stops, their start and end stages, the ticket rates (concession allotted if needed), ticket timing, total count of passengers and the revenue collected in each trip.

4.2 Modelling of Prior Distribution

The steps in modelling the data as prior distribution is shown in Fig. 2.

Figure 2
figure 2

Steps in modelling the prior.

Full size image

The passenger arrival is modelled as a Poisson process as given in Eq. 1. The \(\lambda _i\) values in the Eq. 1 are calculated, based on the minimization of mean square error (MMSE) criterion between the actual statistics and the desired one. If the actual passenger data P(x) and the desired is \(P_d(x)\), the mean square error between P(x) and \(P_d(x)\) is.

$$\begin{aligned} \nonumber \zeta&=E[|P_d(x)-P(x)|]^2 \\&=tr\{E\{(P_d-P)(P_d-P)^T\}\} \end{aligned}$$
(6)

The values of the parameters that minimizes the error is obtained by

$$\begin{aligned} \frac{\partial \zeta }{\partial \lambda _1}=0, \frac{\partial \zeta }{\partial \lambda _2}=0, \frac{\partial \zeta }{\partial \lambda _3}=0, \cdots , \frac{\partial \zeta }{\partial \lambda _n}=0, \end{aligned}$$
(7)

The optimization using Powell algorithm [38], yielded 4 significant \(\lambda\) values. After calculating these values, the arrival process is modelled as a sum of Poisson processes. Then the passenger data is modelled to compute the prior distribution.

4.3 Gaussian Process Modelling

The advantage of Gaussian process regression technique is that it fits best for those data sets with no definite pattern. The traffic data being random in nature, Gaussian process prediction method is the most suitable one. A Gaussian Process is completely specified by a mean and covariance function. The Gaussian Processes is used to perform regression via supervised learning, called Gaussian Process Regression. This regression process is modelled as kernelized Bayesian linear regression, where the kernel parameterization is determined by the choice of covariance or kernel function, as well as the passenger data used to make predictions. Given a training dataset consisting of N observations:

$$\begin{aligned} \textit{D}_{train}=(X,y)=\{x_i,y_i\}_{i=1}^N,x_i,y\in \mathbf {R} \end{aligned}$$
(8)

As well as a test dataset consisting of N’ points:

$$\begin{aligned} \textit{D}_{test}=(X_*)=\{x_{*,i}\}_{i=1}^{N'},x_{*,i}\in \mathbf {R} \end{aligned}$$
(9)

GPR predicts a posterior Gaussian distribution for targets over test points \(X_*\) by computing the parameters of this Gaussian distribution given observed training data. This GPR model is validated and compared with the existing technologies like Student-t process and KRR and the final results are analyzed. The experimental steps for implementing the above methodology are detailed below.

5 Experiment

The experimental steps in passenger prediction along with the tools in performing it are presented in Fig. 3. These steps are implemented using different python modules. The data is read using python pandas module, an open source easy to use tool that is used for handling large data sets. The data is converted as pandas dataframes, which are easily manipulated and cleaned using numpy module in python. The geospatial data is visualized using Basemap, pandas and matplotlib modules. For passenger profiling, tripsheet and passenger reports are read using pandas, analyzed using numpy and plotted using matplotlib module. The prior distribution is computed using numpy, scipy and lmfit modules. After the computation of the prior, GPR is designed and implemented in PyMC3 environment with the support of arviz and seaborn modules. PyMC3, an open source probabilistic programming framework with Theano as backend, allows flexible implementation of Bayesian statistical models in python. The performance of the model is validated using the modules in PyMC3, scipy and sklearn modules.

Figure 3
figure 3

Experimental steps with tools.

Full size image

5.1 Data Analysis

A total of 47 bus routes in Pune city are taken for the study. The complete data related to these set of routes is taken into consideration. The detailed steps in data analysis are shown in Fig. 4. For proper modelling, analysis is categorized into Geospatial analysis and Passenger analysis.

5.1.1 Geospatial Analysis

The geospatial analysis aims in concentrating the land use pattern of the selected routes. The geospatial data of 47 routes include the detailed report of different stops in all routes with their latitude and longitude values. It is plotted to study the geographical extent of the existing routes (Fig. 5a). The latitude and longitude values of each route are read using pandas and converted into individual dataframes. These dataframes are read by the Basemap package to visualize the extent and coverage of individual routes. These routes are plotted in a canvas with the longitude and latitude values of Pune city along horizontal and vertical directions. The visualization is achieved using Basemap and matplotlib module. The plot reveals that these routes cover almost the entire Pune city. To figure out the various statistics and to study the feasibility of different routes, the routes with more than \(80\%\) overlap are identified (Fig. 5b). The identified routes are further studied on the basis of the number of schedules, trip timings and passenger traffic which is detailed in the next section.

Figure 4
figure 4

Steps in data analysis.

Full size image
Figure 5
figure 5

Geospatial visualization of routes.

Full size image

5.1.2 Passenger Analysis

The tripsheet report and passenger profile reports are used for passenger modelling. The steps involved in cleaning and analysis are illustrated in Fig. 6. The data from the ticketing machines needed cleaning since it contained non ASCII characters, missing entries and misplaced values. Approximately more than 500 files are corrected. Python modules are developed for correcting them iteratively. The files of each route are read as pandas dataframes, divided into uptrips and downtrip dataframes using the data obtained from the Tripsheet report. In cases where trips are wrongly numbered as up and down, stages are extrapolated from the passenger profile report to seperate the dataframes correctly into uptrips and downtrips. After separating them, the number of passengers travelling in different stages are calculated for detailed analysis.

Figure 6
figure 6

Steps in passenger traffic analysis.

Full size image

The passenger data analysis is broadly divided into timewise analysis and stagewise analysis of passenger data. The timewise analysis is performed to interpret the travel behaviour of passengers. The travel pattern of passengers vary depending on the day as well as the time of travel. The detailed analysis gives inferences that it shows peak during office hours (Fig. 7a). The day wise analysis (Fig. 7b) is performed to understand passenger flow behaviour and came to a conclusion that it shows peak on working days.

Figure 7
figure 7

Passenger analysis.

Full size image
Figure 8
figure 8

Stagewise passenger data.

Full size image

In order to draw more inferences, the stagewise analysis of passengers is also studied. The Fig. 8 shows the passenger history in route 24, one among the longest and busiest traffic route in Pune city.There are twelve fare stages in this bus route and the graph shows the average number of passengers retained in the bus as the trip proceeds through the different fare stages. It reveals that schedules are running with more than fifteen passengers in almost all the stages and the bus utilization is good. All these analyses leads to the requirement of modelling the passenger arrival.

5.2 Prior Computation

As the arrival of passengers follows a Poisson process, the prior distribution selected is a Gamma Distribution [39]. The Gamma distribution, expressed as

$$\begin{aligned} f(x|\alpha ,\beta )=\frac{\beta ^\alpha }{\varGamma (\alpha )}x^{\alpha -1}e^{-\beta x} \end{aligned}$$
(10)

where \(\varGamma (\alpha )\) is given by

$$\begin{aligned} \varGamma (\alpha )=\int _{0}^{\infty }x^{\alpha -1}e^{-x}dx \end{aligned}$$
(11)

The optimization using Powell algorithm calculated four significant \(\lambda\) values. The gamma function parameters \(\alpha\) and \(\beta\) are calculated using the Poisson rate \(\lambda\) and variance of the data \(\sigma ^2_x\) as

$$\begin{aligned} \alpha&=\frac{\lambda ^2}{\sigma ^2_x} \end{aligned}$$
(12)
$$\begin{aligned} \beta&=\frac{\sigma ^2_x}{\lambda } \end{aligned}$$
(13)

Thus these statistics arrives in choosing Gamma distribution as the prior model in GPR.

5.3 Implementation of GPR

The GPR model is implemented by using PyMC3 environment as it is very flexible and Gamma distribution is the default prior in this environment for GPs. In this environment, GPs are designed and implemented by choosing the covariance function and then calculating the marginal likelihood for prediction. The covariance matrices are often described using functions known as kernels. The one used here is the Exponential quadratic kernel as it results in a smooth prior on functions sampled from Gaussian process. For an input function f(x), it is represented as

$$\begin{aligned} K(x,x')=\eta ^2\exp (-\frac{||x-x'||^2}{2l^2}) \end{aligned}$$
(14)

Here \(||x-x'||^2\) is the squared Euclidean distance given as

$$\begin{aligned} ||x-x'||^2=(x_1-x_{1}')^2+(x_2-x_{2}')^2+\cdots +(x_n-x_{n}')^2 \end{aligned}$$
(15)

\(\eta\) and l represents the hyperparameters to interpret and analyze data. The term l is called the characteristic length scale (or bandwidth or variance) and controls the width of the kernel along a particular axis and \(\eta\) represents the length along other axes or the amplitude. These hyperparameters represent the relevance of each value in the input data. If the l value is large, the covariance becomes almost independent of that particular input data. Hence this data is invalidated and removed, thus implements Automatic Relevance Determination (ARD) process [40]. In our prediction problem, l is specified using Gamma distribution and \(\eta\) using Half normal distribution. The GP is instantiated using this covariance function and marginal likelihood is computed to predict the passenger data. The mean of this distribution is taken as the prediction. To evaluate the model performance, traceplot, posterior plot and autocorrelation of the covariance function parameters are examined.

The route 24 that connects Katraj to Maharashtra Housing Yerwada, covering a distance of \(16\,km\) is selected due to its heavy traffic and it represents the reliable samples of Pune population since it connects the major parts of Pune city. The passenger traffic data for a period of three months (October, November and December 2019) is taken.

5.4 Performance Validation

The performance of GP model is validated using

  • RMSE

  • Kernel density estimation

  • Energy level distribution

The RMSE for the predicted data \(f(x_*)\) is given by

$$\begin{aligned} RMSE=\sqrt{\frac{1}{n}\sum _{n}[f(x_*)-f(x)]^2} \end{aligned}$$
(16)

It indicates the numerical reliability of the model in predicting the traffic data. The kernel density estimation plot is used to compare the pdf of predicted one with the actual one. It is used in validating the predicted passenger data models for choosing proper bus scheduling algorithms. In the energy distribution graph, the difference in energy between the start and the end of the trajectory indicates the efficiency of the sampler. If the distribution shows longer tails, the efficiency of the sampler will deteriorate quickly. The energy levels of the samples are used to identify posteriors with problematically long tails.

The Gaussian Regression model is compared with the existing elliptical processes [41] such as Student-t processes and Kernel Ridge Regression process, predicting the random traffic data in terms of the above parameters. The Student-t has been used in variety of applications [42,43,44,45]. It is derived from the Gaussian process by placing an inverse Wishart process prior on the kernel. The KRR process is used in predicting wind speed [46], molecular orbital energy predictions [47] and in many other applications.

6 Results

The significant \(\lambda\) values calculated, using Powell optimization algorithm are presented in Table 2.

Table 2 Parameters calculated using Powell method.
Full size table

Using these parameters, the passenger data is modelled (Fig. 9) and obtained a Poisson response. The Gamma parameters, computed using these Poisson parameters yielded 10.67 as the \(\alpha\) value and 6.93 as the \(\beta\) value for designing the GPR model.

Figure 9
figure 9

Passenger statistics.

Full size image

The GPR model is used to predict the passenger travel pattern for fifteen days. The gaussian model prediction is shown in Fig. 10a and that by Student-t model in Fig. 10b and by KRR model in Fig. 10c. Unlike the case of Student-t distribution model, the predicted value is almost following the actual one in the case of GPR.

Figure 10
figure 10

Comparison of prediction.

Full size image

The trace plot of the gaussian process model is given by (Fig. 11). The posterior density of noise decreases monotonically. The posterior distribution spread is indicated using Highest Posterior Density (HPD) interval or the bayesian credible interval. It is the shortest interval indicating a given portion of the probability density. The posterior plot, that indicates the distribution of hyperparameter,\(\eta\) and length function, l is shown in Fig. 12. It indicates mean value 45 and HPD of \(94\%\) for the hyperparameter \(\eta\).

Figure 11
figure 11

Traceplot of gaussian model.

Full size image

The precision of the model is benchmarked by a sharply decaying autocorrelation function. An ideal sample from the posterior distribution has autocorrelation value similar to an impulse function. Figure 13 shows that the samples from posterior distribution shows rapidly decaying autocorrelation values, validating the model.

Figure 12
figure 12

Posterior plots.

Full size image
Figure 13
figure 13

Autocorrelation plot.

Full size image

The performance of the model is analyzed, validated and compared in terms of its RMSE value, kde and sampler efficiency.

6.1 Computation of Root Mean Square Error

The mean square error of the predicted data for the GPR model is found to be 5.33, while that for the existing model like Student-t, the RMSE value is 13.36 and for KRR it is 24.61. The performance of the model is much better as the error is decreased by a factor of 8 compared to Student-t and by 19.28 while comparing with KRR model, validating the suitability of the GPR model for predicting passenger traffic.

6.2 Kernel Density Estimation

The statistical properties are analyzed by plotting the joint probability densities (samples and the predicted values) for both models (Fig. 14). In the case of student-t and KRR, it is visible from the figure that the probability density plots of actual and predicted ones are different (Fig. 14b, c) while that with the GPR model (Fig. 14a), both the samples and predicted one follows similar distribution.

Figure 14
figure 14

Kernel density estimation.

Full size image

6.3 Sampler Efficiency

In order to assess the efficiency of the sampler, overall energy level distribution of the samples, with the change of energy between successive samples, are plotted. The sampler efficiency in the case of Student-t is low, as the distribution has longer tails. From Fig. 15, it is clear that GPR sampler (Fig. 15a) is highly efficient than Student-t sampler (Fig. 15b).

Figure 15
figure 15

Energy of the samples.

Full size image

7 Inferences and Conclusion

Gaussian process regression model is designed to predict the bus passenger traffic which shows random nature and high parametric dependence. The GPR model is observed to perform well in terms of prediction accuracy even with limited data. The arrival of passengers at different bus stops follows Poisson distribution whose statistics links with the bus passenger traffic. The passenger data being modelled as a Gamma distribution, along with the likelihood, GPR computes the posterior distribution. Finally the conditional distribution is computed to evaluate the predicted values. The training and prediction process are much simpler and faster than the traditional machine learning techniques. The GPR method took only 170 seconds to complete the training and prediction process compared to 626.46 seconds for the Student-t process. The present work is done with the objective of developing scheduling algorithms. It is also instrumental in designing Intelligent traffic management systems. Such intelligent systems that are adapted to the passenger traffic will improve the financial performance of transportation systems and passenger comfort.