Tailored Wakeby-type distribution for random bus headway adherence ratio

Transportation Research Part C 86 (2018) 220–244

Contents lists available at ScienceDirect

Transportation Research Part C journal homepage: www.elsevier.com/locate/trc

Tailored Wakeby-type distribution for random bus headway adherence ratio

T

⁎

Man Zhanga,b, Qiang Mengb, , Liujiang Kangb, Wenquan Lic a b c

School of Civil Engineering, Suzhou University of Science and Technology, Ke Rui Road #1, Suzhou 215009, China Department of Civil and Environmental Engineering, National University of Singapore, Singapore 117576, Singapore School of Transportation, Southeast University, Si Pai Lou #2, Nanjing 210096, China

AR TI CLE I NF O

AB S T R A CT

Keywords: Bus routeregularity Headway adherence Probability distribution estimation Wakeby distribution L-moment

This paper addresses an interesting and practical bus headway adherence issue for a given public bus route with a number of bus stops. It first defines the random headway adherence ratio (HAR) at a particular bus stop of a specific bus route as the ratio of difference between actual bus headway and scheduled headway with respect to the scheduled headway. This study proceeds to customize a four-step procedure to estimate a probability distribution that can describe the random HAR at each bus stop of the bus route by using the automatic vehicle location (AVL) data. Our real case studies with 44,025 HAR data show that the 19 existing probability distributions including Lognormal, Gamma, Beta and Wakeby are unable to well fit these HAR data. This study thus proposes a tailored Wakeby-type distribution with five parameters. After deriving two fundamental propositions for the tailored Wakeby-type distribution, a tangible L-moment based method to estimate those parameters involved the tailored Wakeby distribution is presented. The tailored Wakeby-type distributions can meet our expectation via our real case studies. Finally, applications of the tailored Wakeby-type distribution derived for the random HAR are conducted.

1. Introduction Public bus as one of the major transport modes plays a vital role in the large Asian cities such as Singapore, Hong Kong, Beijing and Bangkok. Taking Singapore as an example, a city (nation) with a population of 5.3 million people generated around 3.4 million daily bus trips in 2012 that accounted for 50 per cent of the total public transport trips (LTA, 2013). According to the 2012 household interview travel survey of Singapore, about 62 per cent of all commuters’ trips during the peak hour periods were made on public transport. The aim and vision in 2013 Land Transport Master Plane of Singapore are to achieve a peak period public transport mode share of 70 per cent in 2020 and 75 per cent by 2030 by providing more connections and integration between public transport modes - bus, train and other forms of travel like cycling. Nevertheless, Singaporeans' overall satisfaction with the public transport system has dipped from 90.3 per cent in 2011 to 88.8 percent in 2012 and the satisfaction levels for bus and rail reliability also fell, according to the latest public transport customer satisfaction survey conducted in October of 2012. To improve the reliability of public transport, Land Transport Authority (LTA) of Singapore has just implemented a bus service reliability framework (BSRF) for 22 bus routes on a two-year trial since 3 February 2014. Following the quality incentive contract (QIC) of London, the BSRF assesses regularity of bus service by the passengers’ excess waiting time (EWT) that is equal to the actual waiting time (AWT) minus the schedule waiting time (SWT). Incentives and penalties are imposed to the relevant bus operators by

⁎

Corresponding author. E-mail address: [email protected] (Q. Meng).

https://doi.org/10.1016/j.trc.2017.11.013 Received 8 September 2017; Received in revised form 11 November 2017; Accepted 11 November 2017 0968-090X/ © 2017 Elsevier Ltd. All rights reserved.

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

Nomenclature

Hiw,j ,B

Symbols

m n A B [Ts,Te] IA + 1

IB + 1

j

i

τiA τiB hiA hiB w

number of bus stops for the outbound direction of the studied bus route number of bus stops for the inbound direction of the studied bus route departure-terminal for the outbound bus trip departure-terminal for the inbound bus trip time window under concern number of scheduled trips departing from terminal A in time window [Ts,Te] number of scheduled trips departing from terminal B in time window [Ts,Te] bus stop j, j ∈ [1,2,…,m] for outbound direction and j ∈ [m + 1,m + 2,…,m + n] for inbound direction bus trip i, i ∈ [0,1,…,IA] for trips scheduled to depart from terminal A in time window [Ts,Te], i ∈ [0,1,…,IB] for trips scheduled to depart from terminal B in time window [Ts,Te] scheduled departure time of the ith outbound bus trip from terminal A, τiA ∈ [Ts,Te] scheduled departure time of the ith inbound bus trip from terminal B, τiB ∈ [Ts,Te] scheduled bus headways for the ith outbound bus trip from terminal A scheduled bus headways for the ith inbound bus trip from terminal B workday, w ∈ W and

tiw,j ,A tiw,1,A

tiw,j ,B

tiw,m,B+ 1

HARiw,j ,A

HARiw,j ,B Ljw,outbound

Ljw,inbound

λr αr

W= {Monday,Thuesday, Wednesday, Thursday,Friday}

Hiw,j ,A

the outbound bus stop j ( j ∈ {1,2,…,m} ) on a particular workday w ∈ W actual bus headway of the ith trip (i = 1,2,…,IB ) at the inbound bus stop j ( j ∈ {m + 1,m + 2,…,m + n} ) on a particular workday w ∈ W actual arrival time at the outbound bus stop j ( j = 2,3,…,m ) for the ith trip (i = 1,2,…,IA ) on a particular workday w ∈ W Actual departure time from the first outbound bus stop for the ith trip (i = 1,2,…,IA ) on a particular workday w ∈ W actual arrival time at the inbound bus stop j ( j = m + 2,m + 3,…,m + n ) for the ith trip (i = 1,2,…,IB ) on a particular workday w ∈ W actual departure time from the first inbound bus stop for the ith trip (i = 1,2,…,IB ) on a particular workday w ∈ W bus headway adherence ratio at the outbound bus stop j ( j = 1,2,…,m ) for the ith (i = 1,2,…,IA ) trip on a particular workday w ∈ W bus headway adherence ratio at the inbound bus stop j ( j = m + 1,2,…,m + n ) for the ith (i = 1,2,…,IB ) trip on a particular workday w ∈ W the number of actual bus adherence ratios at the outbound bus stop j ∈ {1,2,…,m} on the workday w∈W the number of actual bus adherence ratios at the inbound bus stop j ∈ {m + 1,2,…,m + n} on the workday w ∈ W L-moment of a random variable (r = 1,2,…) probability weighted movements of a random variable (r = 0,1,…)

actual bus headway of the ith (i = 1,2,…,IA ) trip at

comparing the difference between EWT and a predefined baseline EWT because EWT reflects the additional times passenger faces as a result of irregular bus services (LTA, 2014). The EWT at a particular bus stop of a bus route can be characterized by the regularity of actual bus headway at bus stops with respect to the scheduled bus headway, which is referred to as the bus headway adherence (TRB, 2013). The bus headway adherence is, however, determined by a number of factors that are summarized by Liu and Sinha (2007) as follows: (i) road traffic conditions including traffic congestion levels, traffic composition, day-to-day and within-day variation in travel demand and etc., (ii) bus route characteristics including length of a bus route, location of bus stops, provision of bus lanes, number of intersections on the bus route and etc., (iii) passenger characteristics including passenger volume at bus stops, variability in passenger volume and etc., and (v) bus operational characteristics including scheduling system, bus fleet availability, ticketing system, variability in drivers’ behavior and experiences, and etc. Some of these factors have high uncertainties varying by the time of day and day of week. To ensure a rational baseline EWT and an achievable incentive standard, it is thus of great practical significance to investigate the bus headway adherence issue from a probabilistic perspective. A public bus route usually consists of a number of bus stops including two bus terminals, served by a string of buses dispatched from each of these two bus terminals as per a predetermined timetable. The scheduled headway for the bus route is the difference between the departure times of two consecutive buses departing from the same bus terminal. The scheduled headways are uneven in practical even for the peak hours (Ceder, 2007), which calls for a relative index to quantitatively measure the headway adherence and taking both real headway and scheduled headway into consideration. The headway adherence ratio (HAR) is thus proposed to quantify bus headway adherence at a bus stop, which is defined as the ratio of difference between actual headway at the bus stop and the scheduled headway with respect to the scheduled headway for a specific bus trip. A negative (positive) HAR means a percentage of actual headway being less (larger) than the schedule headway. The closer HAR is to zero the better performance of the bus regularity. Above all, the HAR makes performance of the headway adherence comparable for those bus routes with different scheduled headways. It is more reasonable to assume that the HAR is a random variable, called the random HAR hereafter, due to the uncertainty of the aforementioned factors affecting the bus headway adherence. Such an assumption can be further validated by the stochastic bus dwell time at a bus stop (Meng and Qu, 2013; Bian et al., 2015) and random bus travel time (in-vehicle travel time) between two consecutive bus stops (Polus, 1975, 1979). This study aims to derive a type of probability distribution that can 221

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

appropriately describe the random HAR at each bus stop of a bus route in peak hours by using the AVL data. The derived probability distribution not only enables a better understanding of bus headway regularity, it is also essential to deal with bus headway related research issues including bus bunching control (e.g. Daganzo, 2009; Daganzo and Pilachowski, 2011; Hickman, 2001; Munoz et al., 2013; Yu et al., 2016), passenger waiting time modelling (Chang and Hsu, 2004), transit assignment (e.g., Szeto et al., 2013; Wu et al., 1994; de Cea and Fernández, 1993), bus timetable optimization (e.g. Ceder et al., 2001; Ibarra-Rojas and Rios-Solis, 2012; Yan and Tang, 2008), bus route simulation (Cortes et al., 2010; Toledo and Cats, 2010; Petersen et al., 2013) and coordination between bus and train services (Dou et al., 2015; Kang and Meng, 2017). 1.1. Relevant studies Our throughout literature review via the academic search engines – Science Direct, Scopus, Web of Science and Google Scholar – shows that the random HAR proposed by this study is a new concept and it has not yet been investigated in literature. However, there are two categories of relevant studies: bus headway distribution and bus headway adherence measurement, and they are deserving of review as follows. For the bus headway distribution, a few research works can be found in literature. Polus (1975) is one of the early studies to identify a probability distribution for the random bus headway at a bus stop by using the statistical analysis methods. However, he used the bus travel time data instead of bus headway data. Actual bus headway under this circumstance was defined as the sum of dispatching headway at the bus terminal and difference between travel times of two successive bus arrivals at a bus stop. By assuming one constant for all the dispatching bus headways, a probability density function (PDF) of bus headway can be derived from that of the random bus travel time. The Normal, Gamma and Beta distributions were examined for the bus travel time distribution, and the good of fitness (GoF) test using the Kolmogorov-Simonov (K-S) statistic with the significance level of 0.2 was conducted. According to the test results, the Beta distribution was chosen as the best-fitted probability distribution for the random bus travel time. Hence, the bus headway distribution can be derived from the difference of two independent and identical Beta distributions, which was named as bd-Beta distribution by Polus (1975). Note that the bd-Beta distribution is far away from the conventional Beta distribution. Polus (1979) also validated the bd-Beta distribution using the bus headway data collected from several bus routes in Evanston, Illinois, a Chicago suburb, and concluded that the bd-Beta distribution with the calibrated parameters is able to formulate the random bus headway. Adebisi (1986) explicitly defined that the random bus headway at a specific bus stop of a bus route comprises the random bus travel time from the preceding bus stop to the bus top and random bus dwell time occurred in the preceding bus stop. He proposed an explicit formula to estimate variance of random bus headway for each bus stop based on the three unrealistic assumptions: (i) all the dispatching headways are the same (ii) the unit bus travel time between any two consecutive bus stops follows the same probability distribution (iii) passengers arriving at a bus stop follows a stationary Poisson process. Computer simulations have been also used to examine some characteristics of bus headway at a bus stop (Estrada et al., 2016; Wu et al., 2017). Nagatani (2001) proposed a bus headway simulation model by assuming that (i) “a bus driver operates his bus in such a manner that his velocity increases or decreases according as the headway is large or small” and (ii) “the mean velocity of a bus j at bus stop m depends only on the headway of bus j at bus stop m”. These two assumptions are obviously problematic from public transit engineering point of view. As the proposed bus headway distribution is a simple application of the car-following models in traffic flow theory, stability analysis for the bus headway is conducted and some naive conclusions are drawn. Such a research methodology has been also adopted by Hill (2003) and Nagatani (2003) for the bus headway studies. Of course, these studies can only exhibit dynamics of bus headway rather than its probability distribution. Bellei and Gkoumas (2010) developed a stochastic simulation model to examine the bus headway distribution by mainly assuming that (i) the bus travel time between two adjacent bus stops followed a triangular distribution; (ii) alighting and arriving passengers at a bus stop followed Binomial and Poisson distribution, respectively; (iii) bus dwell time was merely proportional boarding and alighting. To implement such a simulation model, numerical values of these assumed distributions were specified. With 1000 simulations cycles, headway histograms of the main stops were obtained. Aiding with the curve fitting toolbox of Matlab, based on the method of Generalized Least Squares, the representative points, which have the intermediate point of each headway length interval as abscissa and the percentage of headways whose length is within the interval as ordinate, were interpolated using first and second degree Gauss curves. Finally, the adjusted R-square values indicated that the first degree Gauss curves are sufficient to depict the situation at initial stops, while the second Gauss curves are needed to depict the headway distribution at the final stops. Toledo and Cats (2010) put forward a novel mesoscopic simulation model for transit network operations and they did not explore the bus headway distribution. The widespread application of Automatic Vehicle Location (AVL) technology in public transport systems has opened a new way for estimating a probability distribution of the random bus headway. Strathman et al. (1999) realized randomness of the bus headway after analyzing the AVL data from several bus routes in Portland, Oregon. However, they could not analyze the probability distribution of the random bus headway. Lin and Ruan (2009) defined an interesting probability-based headway regularity measure using probability distribution of bus dispatching headway. They used one-month AVL data in the morning peak hours (7:00 am–9:00 am) from one bus route in Chicago to estimate the probability distribution of bus dispatching headway. Gamma, Weibull, Logistic and Lognormal distributions were fitted on the collected dispatching headways and the goodness-of-fit statistics of K-S, Anderson-Darling (A-D) and Chi-square showed that the dispatching headway of the studied bus route follows the Gamma distribution. This finding also demonstrates that the constant dispatching headway assumed in the aforementioned studies is unrealistic. Nevertheless, histograms of the bus headways at the other locations shown by this study indicate that the calibrated Gamma distribution is unfit for some locations. Based on the AVL data, Guo et al. (2011) and Saberi et al. (2013) analyzed the bus headway using the percentage rank rather than the bus headway distribution identification. For more data-driven techniques, interested readers refer to Byon et al. 222

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

(2017) and Nesheli and Ceder (2016). Without resorting to the bus headway data, it is almost impossible to derive an analytical probability distribution of the random bus headway under some mild and realistic assumptions due to those uncertain factors involved in the bus operations. Compared to the computer simulation models, identifying the bus headway distribution using the AVL data is more appropriate although such an approach has not well documented in literature. As for the bus headway adherence measurement, several indexes or metrics have been proposed or used to measure bus headway adherence for a particular bus stop of a specific bus route, out of which three representatives are illustrated below. The first one called the percentage regularity deviation mean (PRDM) is defined by averaging the absolute values of ratio of difference between the actual and scheduled headway to the scheduled headway for those bus trips occurred in an interested time period, as shown below: n

∑i = 1 PRDMj =

Actual Headwayi,j − ScheduledHeadwayi,j Scheduled Headwayi,j

(1)

n the ith

actual and where j stands for the target bus stop, n is the number of bus trips, Actual Headwayi,j and Scheduled Headwayi,j denote scheduled bus headways at bus stop j , respectively, where i = 1,2,…,n . The PRDM is originally proposed by Hakkesteegt and Muller (1981), and it has been subsequently utilized by Ap. Sorratini et al. (2008), van Oort and van Nes (2009) and van Oort (2011) for transit network reliability analysis. The second one is the percent of schedule deviation (PSD) defined by actual departure minus scheduled departure that fall within a defined range to measure the on-time transit service performance for a bus stop (TRB, 2013). The third one is the coefficient of variation (CoV) of headway defined by the standard deviation of headway divided by the mean headway for a particular bus stop, proposed by Adebisi (1986) and compiled by TRB (2013). It can be seen that both PRDM and PSD are the statistical metrics calculated from sample headway data. Hence, they lack of capability to probabilistically analyze the bus headway adherence issue. The CoV of headway acknowledges randomness of bus headway although it is estimated from headway sample data without resorting to the probability distribution of the random bus headway. 1.2. Objectives and contributions There is no doubt that a probability distribution of the random HAR will provide an even richer way for the bus reliability analysis than any bus headway adherence metrics estimated from the HAR data. Although the random HAR is highly correlated to the bus headway distribution, its probability distribution is most likely unable to directly derive from the bus headway distribution because of the uneven scheduled bus headways at a bus terminal in practice. To provide a comprehensive understanding of bus headway performance and lay a foundation for bus service reliability analysis, the objective of this study is to estimate a type of probability distribution that can well describe the random bus HAR at most bus stops of a particular bus route during the peak hours by using the AVL data. As in much research such as bus operation simulation or transit assignment, distribution about bus headway needs assumptions, and only which fits most bus stops is valid and efficient. For the sake of presentation, such a type of probability distribution is called a unified probability distribution. To do so, we will first customize a four-step procedure to estimate a unified probability distribution using the AVL data. The customized four-step procedure is applied for three bus routes in Changzhou city of China. It will be found that the existing probability distributions are unable to well fit the bus HAR data of the three bus routes. We thus proceed to derive a novel probability distribution, referred to as the tailored Wakeby-type distribution, which fits the collected bus HAR data much better than the existing probability distributions. Finally, application of the tailored Wakeby-type distribution derived for the random HAR will be discussed. The contributions of this study with respect to the current literature and industrial applications can be summarized as follows: (a) The random HAR defined by this study is a new concept with significant practice relevance. The defined random HAR is a novel concept according to our literature review and its probability distribution can provide better understanding of bus headway adherence. More importantly, it can be used to determine an appropriate bus headway control strategy or evaluate effectiveness of some bus bunching control strategies. It can also be used to estimate a base EWT in the Bus Service Reliability Framework (BSRF) implemented by the LTA of Singapore. In reality, this concept was initially created when we had a technical meeting on determination of the base EWT with a traffic engineer from the LTA of Singapore. (b) Usefulness of the customized four-step procedure to identify a unified probability distribution of the random HAR. AVL data are nowadays available in the public transport industry. The customized four-step procedure provides a useful and tangible tool for traffic engineers to estimate a probability distribution of the random HAR by using the AVL data. This procedure integrates a few statistical analysis and parameter estimation techniques. (c) Novelty of the tailored Wakeby-type distribution. The proposed tailored Wakeby-type distribution is a new probability distribution and it is creatively developed based on the Wakeby distribution. It can play the unified probability role for the random HAR in accordance with our case studies. A useful L-moments based method to estimate those parameters involved in the tailored Wakeby-type distribution is also presented. The remained of this paper is organized as follows. Section 2 gives notations and assumptions used by this study, and the problem description. Section 3 presents a customized four-step procedure to estimate a unified probability distribution for the random HAR. Section 4 proposes a tailored Wakeby-type distribution with five parameters followed by two fundamental propositions a tangible Lmoment based method to estimate those parameters involved the tailored Wakeby distribution. Section 5 illustrates applications of 223

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

the tailored Wakeby-type distribution of the random HAR. Conclusions are finally presented in Section 6. For readability, some of notations used by this study are presented in Appendix A. 2. Notation, assumptions and problem description The itinerary of a typical public bus route forms a directed loop in terms of bus stops visited by a bus, including two bus terminals and two travel directions differentiating by these two terminals: outbound and inbound. For the sake of presentation, these two bus terminals are named by A and B. A public bus route with m + n bus stops can be schematically represented by its bus stop visiting outbound inbound      sequence, 1 → 2 → ⋯→m → m + 1 → m + 2 → ⋯→m + n → 1, where the numbers 1 and m denote the first and last bus stops of call by a bus for the outbound direction; the numbers m + 1 and m + n denote the first and last bus stops of call by a bus for the inbound direction; bus stops 1 and m + n are in terminal A; and bus stops m and m + 1 are in terminal B, as shown in Fig. 1. The bus route shown in Fig. 1 is served by a string of buses following a daily departure time schedule at each bus terminal. Without loss of generality, it is assumed that these daily scheduled departure times are the same for each of the five workdays grouped into the set W= { Monday, Thuesday, Wednesday, Thursday, Friday} . Since this study focuses on the bus headway adherence analysis in the morning or afternoon peak hours, we choose a time window denoted by [Ts,Te] over which bus trips scheduled to dispatch from both terminals can capture the morning or afternoon perk hours. The daily scheduled bus departure times covered by the time window can be sorted from the earliest one to the latest one for each of the two terminals, as shown in Table 1. In Table 1, τiA (i = 0,1,…,IA ) denotes the departure time of the ith outbound bus trip from terminal A and τiB (i = 0,1,…,IB ) denotes the departure time of the ith inbound bus trip from terminal B. It should be noted that the total numbers of bus trips from both terminals may not be the same, namely, (IA + 1) may not be equal to (IB + 1) . Table 1 transparently shows that there are IA and IB scheduled bus headways, denoted by hiA (i = 1,2,…,IA ) and hiB (i = 1,2,…,IB ), for the outbound and inbound directions, respectively, which can be calculated by:

hiA = τiA−τiA− 1,i = 1,2,…,IA, and

(2)

hiB

(3)

=

τiB−τiB− 1,i

= 1,2,…,IB.

Note that these scheduled bus headways may not be the same in practice, i.e., uneven scheduled bus headways. Since the road traffic conditions may vary from one workday to another one, the actual bus headway at a particular bus stop of the bus route, which is defined by the difference of arrival times at the bus stop by two consecutive buses serving the bus route, should be differentiated by each workday. Let Hiw,j ,A be the ith (i = 1,2,…,IA ) actual bus headway at the outbound bus stop j ∈ {1,2,…,m} on a particular workday w ∈ W , and Hiw,j ,B be the ith (i = 1,2,…,IB ) actual bus headway at the inbound bus stop j ∈ {m + 1,m + 2,…,m + n} on a particular workday w ∈ W . These actual bus headways can be calculated by A Hiw,1,A = tiw,1,A−tiw−,1,1 ,i = 1,2,…,IA;w ∈ W

(4)

Hiw,j ,A = tiw,j ,A−tiw−,1,Aj,i = 1,2,…,IA;j = 2,…,m;w ∈ W

(5)

Hiw,m,B+ 1 = tiw,m,B+ 1−tiw−,1,B m + 1,i = 1,2,…,IB;w ∈ W

(6)

Hiw,j ,B

(7)

=

tiw,j ,B−tiw−,1,B j,i

= 1,2,…,IB;j = m + 2,m + 3,…,m + n;w ∈ W

where

tiw,j ,A: actual arrival time at the outbound bus stop j exclusive of the first outbound bus stop, namely, j = 2,3,…,m , by the bus dispatched for the ith trip (i = 1,2,…,IA ) from terminal A on a particular workday w ∈ W ; tiw,1,A: actual departure time from the first outbound bus stop (i.e., terminal A) by the bus dispatched for the ith trip (i = 1,2,…,IA ) from terminal A on a particular workday w ∈ W ; tiw,j ,B: actual arrival time at the inbound bus stop j exclusive of the first outbound bus stop, namely, j = m + 2,m + 3,…,m + n , by the bus dispatched for the ith trip (i = 1,2,…,IB ) from terminal B on a particular workday w ∈ W ;

Fig. 1. A Schematic Representation of the Public Bus Route.

224

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

Table 1 Daily Scheduled Bus Departure Times in Time Window [Ts,Te]. Terminal A

Terminal B

Bus trip no.

Departure time

Bus trip no.

Departure time

0

τ1A

0

τ1B

1

τ2A ⋮

1

τ2B ⋮

⋮ IA

⋮ IB

τIAA

τIBB

tiw,m,B+ 1: actual departure time from the first inbound bus stop (i.e., terminal B) by the bus dispatched for the ith trip (i = 1,2,…,IB ) from terminal B on a particular workday w ∈ W .Based on the scheduled bus headways shown in Eqs. (2) and (3) as well as the actual bus headways expressed by Eqs. (4)–(7), we define the bus headway adherence ratio (HAR) for a given bus trip at a specific bus stop of the bus route on a particular workday as follows: HARiw,j ,A =

HARiw,j ,B =

Hiw,j ,A−hiA hiA Hiw,j ,B−hiB hiB

,i = 1,2,…,IA;j = 1,2,…,m;w ∈ W

(8)

,i = 1,2,…,IA;j = m + 1,2,…,m + n;w ∈ W

(9)

The defined bus headway adherence ratio is able to quantify the bus headway adherence issue at a bus stop of a bus route because it conveys the percentage deviation with respect to the scheduled headway. For example, a negative value of HAR implies the earlier arrival; the closer the HAR is to 0, the higher level is the headway adherence. Also, the non-negativity of the bus headway implies that

HARiw,j ,A ⩾ −1,i = 1,2,…,IA;j = 1,2,…,m;w ∈ W

(10)

HARiw,j ,B ⩾ −1,i = 1,2,…,IB;j = m + 1,2,…,m + n;w ∈ W

(11)

However, in reality, random disturbance due to road traffic conditions including congestion and traffic signals, and inhomogeneous driving behavior of bus drivers cause variance to the bus headway adherence. It is thus more rational to formulate the bus HAR at a specific stop j ( j = 1,2,…,m + n ) on a particular workday w ∈ W as a random variable, denoted by RHARjw . As a random

Fig. 2. Four-step Procedure for the Unified Probability Distribution Estimation.

225

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

variable, RHARjw is defined by its probability space (Ω,F jw,Pjw ) with the sample space Ω = [−1,+∞) due to Eqs. (10) and (11), σalgebra F jw and the probability measure Pjw . In other words, RHARjw can be described by its probability density function (PDF) or cumulative density function (CDF) defined on the sample space Ω = [−1,+∞) . This study aims to find a suitable probability distribution that is able to describe the random headway adherence ratio associated the bus trips dispatched during the given window [Ts,Te] for each bus stop of the bus route shown in Fig. 1 on each workday. To achieve this objective, the HAR sample data related to each bus stop are required, denoted by {HARiw,j ,A,l : l = 1,2,…,Ljw,outbound } , where Ljw,outbound is the number of actual bus adherence ratios at the outbound bus stop j ∈ {1,2,…,m} on the workday w ∈ W , or {HARiw,j ,B,l : l = 1,2,…,Ljw,inbound } , where Ljw,inbound is the number of actual bus adherence ratios at the inbound bus stop j ∈ {m + 1,2,…,m + n} on the workday w ∈ W . Based on these HAR sample data, we can estimate a probability distribution of the

Table 2 19 Candidate Probability Distributions. Distribution

Probability density function (PDF)

Parameters to be estimated

Est.M

Sample space Truncated Normal

f (x |μ,σ ,ξ ) =

x−μ 1 ⎞ ϕ⎛ σ ⎝ σ ⎠ ξ − μ⎞ ⎛ 1−Φ ⎝ σ ⎠

f (x |α,β,ξ ) =

α β

⎜

Log-Logistic Erlang Inverse Gaussian Fatigue Life

f (x |β,m,ξ ) =

⎛1 + ⎝

x−ξ β

exp ⎛− ⎝

x−ξ β

β x−ξ

−

Weibull Frechet Pert

f (x |α,β,ξ ) =

⎠

×

⎟

exp(−β / (x − ξ ))

α β

f (x |α,β,ξ ) =

α β

( ) ( ) x−ξ β β x−ξ

α−1

α+1

x−ξ β

β x−ξ

=

f (x |α,β,ξ ) =

f (x |α,β,ξ ) =

α−1 x − ξ⎞ αk ⎛⎜ ⎟ ⎝ β ⎠ α k+1 x − ξ⎞ ⎞ ⎛ β ⎜1 + ⎛⎜ ⎟ ⎟ β ⎝ ⎠ ⎠ ⎝

Johnson SB

f (x |δ ,λ,γ ,ξ ) =

MLM

MLM

α,β

MLM

α

ξ
MLM MLM MLM

5a − ξ − 4m a−ξ

(x − ξ )α − 1 exp(−(x −ξ )/ β ) β α Γ(α )

Burr

α

( ) ⎞⎠ exp ⎛−( ) ⎞⎠ ⎝ exp ⎛− ⎝

4m + a − 5ξ ,α2 a−ξ

MLM

μ,σ ξ
ξ
(x − ξ )α1− 1 (a − x )α2− 1 1 B (α1,α2) (a − ξ )α1+ α2− 1

f (x |m,a,ξ ) =

MLM

ξ
β Γ(α )((x − ξ ) / β )α + 1

where α1 = Gamma

β,m (Integer) ξ⩽x<+∞ |λ,μ ξ⩽x<+∞ α,β

2 ⎛ 1 ln(x − ξ ) − μ ⎞ ⎞ exp ⎜− ⎛ ⎟ 2⎝ σ ⎠ ⎠ ⎝ (x − ξ ) σ 2π

f (x |α,β,ξ ) =

MLM

)⎞ ⎠

⎜

f (x |μ,σ ,ξ ) =

α,β ξ⩽x<+∞

λ (x − ξ − μ)2 ⎞ 2μ2 (x − ξ )

2α (x − ξ )

1

Pearson Type 5

α −2

( ) ⎞⎠

(x − ξ ) / β + β / (x − ξ )

ϕ⎛ ( ⎝α Lognormal

( )

α−1

λ 2π (x − ξ )3

f (x |α,β,ξ ) =

MLM

⎟

x−ξ β

(x − ξ )m − 1 exp(−(x −ξ )/ β ) βmΓ(m)

f (x |λ,μ,ξ ) =

μ,σ ξ⩽x<+∞

1 δ exp ⎜⎛− 2 λ 2π z (1 − z )

⎝

α,β ξ⩽x<+∞ α,β

MLM MLM

ξ⩽x<+∞

2

(γ + δln ( ) ) ⎞⎠ z 1−z

⎟

x−ξ λ kα − 1 k (x − ξ ) exp(−((x −ξ )/ β )k ) βkα Γ(α )

δ ,λ,γ ξ⩽x⩽ξ+λ

MoM

α,β,k ξ⩽x< a,α1,α2 ξ⩽x⩽ α1,α2,β ξ⩽x< a,α1,α2 ξ⩽x< λ1,λ2,γ ξ⩽x<

MLM

where z ≡

Generalized Gamma Beta Pearson Type 6 Kumaraswamy

f (x |α,β,k ,ξ ) =

f (x |a,α1,α2,ξ ) =

(x − ξ )α1− 1 (a − x )α2− 1 1 B (α1,α2) (a − ξ )α1+ α2− 1

f (x |α1,α2,β,ξ ) =

((x − ξ ) / β )α1− 1 βB (α1,α2)(1 + (x − ξ ) / β )α1+ α2

f (x |a,α1,α2,ξ ) =

α1 α2 z α1− 1 ( 1 − z α1 )α2− 1 Where (a − ξ )

Phased Bi-Exponential

f (x |λ1,λ2,γ ,ξ ) = Wakeby

z≡

x−ξ a−ξ

λ1 e−λ1 (x − ξ ) ξ⩽x⩽γ ⎧ ⎨ λ2 e−λ2 (x − γ ) − λ1 (γ − ξ ) γ < x < + ∞ ⎩

x [F |α,β,δ ,γ ,ξ ] = ξ +

γ α [1−(1−F ) β]− [1−(1−F )−δ ] β δ

where F ∈ [0,1)

226

+∞ MLM

a MLM

+∞ MLM

+∞ LSM

+∞

α,β,δ ,γ ξ ≤ x ≤ ∞ (if δ ≥ 0) ξ ≤ x ≤ ξ + α/β − γ/δ (if δ < 0 or γ = 0)

L-M

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

random HAR using the probability and statistical analysis methods. The estimated probability distribution can allow the operators to analyze how often the headway adherence ratios fall in a desired threshold from a probabilistic point of view. For instance, the Massachusetts Bay Transportation Authority (MBTA) uses the percentage of headway greater than 1.5 scheduled headway as an indicator of bus service quality (TRB, 2006). A suitable probability distribution for the random HAR also aids making well-grounded decisions on the bus headway performance measure or proper bus control strategies. 3. Probability distribution estimation for random HAR It is an effective way to estimate a unified probability distribution that can well describe the random HAR at each bus stop of a bus route by means of the probability and statistical analysis methods based on the HAR data: {HARiw,j ,A,l : l = 1,2,…,Ljw,outbound } and {HARiw,j ,B,l : l = 1,2,…,Ljw,inbound } , where j ∈ {1,2,…,m + n} . These data can be obtained from the AVL system installed on buses. Based on several tangible data processing methods (e.g. Mcleod, 2007), the sample date with high effectiveness and accuracy could be ensured. Thus, we first customize a four-step procedure elaborated below to estimate a unified probability distribution using these data. 3.1. Four-step procedure to estimate the unified probability distribution Fig. 2 illustrates the customized four-step procedure using the HAR data to fit a unified probability distribution. The purpose of Step 1 shown in Fig. 2 is to find all the possible candidate probability distributions for the random HAR. This can be done by plotting histogram of the HAR data {HARiw,j ,A,l : l = 1,2,…,Ljw,outbound } and {HARiw,j ,B,l : l = 1,2,…,Ljw,inbound } for each bus stop j ∈ {1,2,…,m + n} . We can subjectively inference some possible probability distributions for further analysis according to these histograms. We have screened all the existing probability distributions and found the 19 possible candidates shown in Table 2 that may be suitable for the random HAR. We should set the location parameter ξ = −1 for all these 19 candidate probability distributions because the sample space of the random HAR is [−1,∞) . The 19th (last) probability distribution listed in Table 2 is called Wakeby distribution without analytical PDF or CDF. It is defined by the quantile function (i.e., inverse cumulative probability function):

x [F ] = ξ +

γ α [1−(1−F ) β]− [1−(1−F )−δ ] ,∀ F ∈ [0,1) δ β

(12)

where parameters α,β,δ ,γ should fulfill the following conditions:

β ≠ 0,δ ≠ 0

(13)

γ⩾0

(14)

β+δ>0

(15)

α+γ⩾0

(16)

Hosking (1986) gave a set of conditions for parameters α,β,δ ,γ to ensure that x (F ) is a quantile function of probability distribution with finite mean, which are presented below: Proposition 1. x (F ) is the quantile function of a probability distribution with finite mean if

β > −1 and β ≠ 0

(17)

δ < 1 and δ ≠ 0

(18)

β+δ⩾0

(19)

α+γ⩾0

(20)

γ⩾0

(21)

Wakeby distribution was introduced by Houghton (1978) for modeling flood flows. It has been successfully used in hydrology to characterize probability distributions of rain fall, flood frequency and so on. Wakeby distribution can describe a wider range of distributional shapes better than those well-recognized probability distributions. With suitable values of parameters α,β,δ ,γ , it is able to mimic shape of many skew distributions such as Extreme-Value, Lognormal and Log-Pearson Type III (Hosking 1986). Interestingly, Wakeby distribution was also applied to examine the bus dwell time by Rashidi and Ranjitkar (2013). Step 2 shown in Fig. 2 aims to estimate parameters involved in a candidate probability distribution by the statistical parameter estimation methods using the HAR data, including the moments based methods (MoM), the maximum likelihood method (MLM) and the least squares method (LSM). Detail of these parameter estimation methods can be found in Ross (2009) and Washington et al. (2011). These conventional parameter estimation methods are, however, inapplicable to estimate the four parameters of Wakeby distribution because Wakeby distribution does not have an analytical PDF or CDF. Hosking (1990) proposed the L-moments based method to estimate these parameters. The appropriate parameter estimation method for each of the 19 probability distribution is listed in the last column of Table 2 because it is suggested by several statistical analysis software packages such as Easyfit (http:// www.mathwave.com/easyfit-distribution-fitting.html) and Bestfit (http://www.palisade.com/devkits/bdk.asp). 227

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

Step 3 shown in Fig. 2 examines how well a candidate probability distribution with estimated parameters fits the HAR data by means of the goodness-of-fit (GoF) test methods. The most common used GoF test statistics include Kolmogorov-Smirinov (K-S), Anderson-Darling (A-D), and Chi-squared ( χ 2 ). Note that the Chi-squared ( χ 2 ) GoF test statistic is appropriate for the large sample size in order to ensure the chi-squared approximation to be valid. As for the A-D test statistic, the critical values for the fitted distribution should be available; however, it is not the case for many candidate distributions shown in Table 2. We therefore suggest using the K-S or Chi-squared ( χ 2 ) test statistics to perform the GoF test with a given significant level denoted by SL. Step 4 shown in Fig. 2 attempts to choose an appropriate probability distribution as the unified probability distribution can well formulate the random HAR for each bus stop j ∈ {1,2,…,m + n} on each workday w ∈ W according to the GoF test results. In principle, at a bus stop level, a probability distribution with the lowest statistic value is the best fitting one for that particular bus stop. However, due to the natural degradation of headway adherence at bus stops lying along the bus route, the most fitted distributions for different bus stops may not be the same. In order to provide a general conclusion to analysts, at a route level, a unified probability distribution is our objective, with proper parameters which could characterize the performance of random HAR in most of the cases. With these two regards, the objective distribution could be preliminarily identified by the fewest rejection times in GoF tests. Moreover, reasonability and applicability are also need to be taken into consideration in the final decision making. It is an advantage if the structure of the distribution model is based on explicit theoretical reasoning about the characteristics of bus operations. The parameters of such model can give additional information on the properties of headway performance. In all, the final decision on the best fitted distribution should take many factors into consideration comprehensively.

3.2. Testing of the Four-step procedure To assess applicability of the four-step procedure to estimate a unified probability distribution suitable for most bus stops of a bus route, we have collected the actual bus headway data at each bus stops of three public bus routes - B1, B11 and B23 (as shown in Fig. 3) – in Changzhou, a middle-sized city with a population of 3.6 million and land mass of 4,375 square kilometers in the eastern China. Public bus services in Changzhou city have been highly developed. By the end of 2012, there were 197 bus routes operating in the urban district, with nearly1.2 million daily riders, which accounted for more than one quarter of the total daily travel trips in Changzhou (Changzhou Statistics, 2013). As shown in Table 3, B1 is a bus rapid transit (BRT) route with advanced BRT stations and dedicated bus lanes. B23 is a traditional

Fig. 3. Route of bus lines B1, B11 and B23.

228

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

bus route where buses are operated in a mixed-traffic environment in the sense that buses and social vehicles travel in the same lane. B11 is a mixed bus route in which only some of its route sections possess dedicated bus lanes, and it can be redeemed as a mix of the BRT route and traditional bus route. Given the above, B1, B11 and B23 could be considered as representatives of three typical bus routes: pure BRT route, traditional bus route, and mixed bus route, which cover almost all the common types of bus routes. In addition, all buses (or vehicles) serving these three routes are equipped with the AVL system. No time-check bus stops are set up along these three routes, and drivers are only expected to endeavor to maintain the scheduled arrival time at the end of each trip. The actual bus headways at each bus stop are collected during the four-week period from 2 September 2013 to 29 September 2013 with bus dispatching time window [6:00AM, 8:00AM], namely Ts = 6: 00AM and Te = 8: 00AM, which fully covers the morning peak hours of Changzhou city. We downloaded the AVL data from the Automated Data Collection (ADC) system of operator of the three bus routes - No.3 Bus Company of Changzhou affiliated to Changzhou Transit Group - for all the buses deployed on these three bus routes after their full-day services. These downloaded AVL data include arrival and departure times of each bus at every bus stop. Data processing is the first step before they can be utilized in our study. Details of the data processing methods can be found from the studies made by Barabino et al. (2013), Ehrlich (2010) and McLeod (2007). With the processed AVL data, the corresponding actual bus headways can thus be calculated according to Eqs. (4)–(7). However, primarily for simplicity and readability of exposition, only the inbound direction shown in Table 3 are selected in the case studies. As a rule of thumb, HAR of different route types might have different probability distributions, even on different workdays. We dealt the data of different bus routes on different workdays separately. Table 4 gives the scheduled bus departure times at the corresponding terminals for the inbound directions of these three bus routes. Based on the collected actual bus headway data and the scheduled bus headway shown in Table 4, the corresponding HAR data at each bus stop on each workday for each of these three bus routes can be obtained according to Eqs. (8) and (9). As a result, we have collected 44,025 HAR data in total, which are summarized in Table 5. Following the four-step procedure shown in Fig. 2, we can estimate a unified probability distribution for the random HAR of bus stops for each bus route on each workday by using the collected HAR data. Let us take the inbound direction with 30 bus stops for the bus route B1 on Monday as an example. To figure out a possible candidate probability distribution, we plot the histograms of the HAR data at each of these 30 bus stops (from the departure stop, numbered 1, to the terminal stop, numbered 30, in sequence) on the Monday, as shown in Fig. 4. Nevertheless, it is very hard for us to “guess” a good candidate probability distribution as per these 30 histograms because of their diversity. It can be also speculated that a unified and well fitted probability distribution is unlikely to be those probability distributions with one or two parameters, such as normal and exponential distributions, due to their limited capability to describe both the sharp peak and long tail exhibited in Fig. 4. To verify this conclusion, we make a trial test for each probability distribution listed in Table 2. We first estimate the parameters involved in each of these 19 probability distributions by the appropriate parameter estimation method shown in the column entitled “Est.M” of Table 2. We adopt the K-S test statistic at a significance level (SL) of 0.05 for the GoF test. The number of rejections for the null hypothesis that the HAR data are drawn from a chosen probability distribution is tabulated in Table 6. Table 6 shows that the number of rejections for these 19 probability distributions is far away from our expectation, which could not be explained by chance. In addition, fitting results of these distributions are random, that is, a rule such as “typeⅠdistribution fits the first stops well and typeⅡdistribution fits the latter stops well” is hard to be found. Fitting results for the other HAR samples of other routes on other workdays are almost the same, and numbers of rejections are also far beyond the expectation. However, to all HAR data of all bus routes on all workdays, Wakeby distribution is outperformed by its fewest rejections, and thus deserves a depth study. The numbers of rejections of Wakeby distribution for all the HAR data from the three bus routes on each workday are shown in Table 7. One interesting finding in Table 7 is that: the number of rejections for route B1 is larger than routes B11 and B23. To explain this result, characteristics of route B1 should be taken into account. B1 is a backbone of the Changzhou BRT network. It stretches from north to south, and covers many high demand distribution points, such as bus terminals, railway stations, CBD, universities, hospitals and major tourist attractions. Seen from Table 3, the length of route B1 is significantly higher than B11 and B23. More potential sources of the uncertainty of the factors affecting the bus headway is resulted by all of these characteristics which are most likely to affect bus headways. Hence, route B1 has the most rejections. 4. Tailored Wakeby-type distributions for the random HAR and their parameter estimation We propose the following tailored Wakeby-type distribution with five parameters for the random HAR with the sample space [ξ0,∞) , where location parameter ξ0 is given, by adding one parametric term in Wakeby distribution expressed by Eq. (12). [Tailored Wakeby-type distribution with five parameters]: Table 3 Characteristics of Three Bus Routes. Route Type

BRT Route Mixed Bus Route Traditional Bus Route

Bus Route No.

B1 B11 B23

Route Length (km)

29.5 18.7 13.7

Number of Bus Stops

Average Operating Speed (km/h)

Inbound Direction

Outbound Direction

30 28 31

30 29 30

229

22.41 19.34 16.44

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

Table 4 Scheduled Headways for the Inbound Direction in Time Window [6:00 am, 8:00 am] B1

B11

B23

Trip No.

Scheduled Departure Time

Scheduled Headway (min)

Trip No.

Scheduled Departure Time

Scheduled Headway (min)

Trip No.

Scheduled Departure Time

Scheduled Headway (min)

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

6:00 6:05 6:09 6:16 6:20 6:28 6:31 6:35 6:42 6:46 6:50 6:53 6:56 6:59 7:04 7:07 7:11 7:14 7:17 7:20 7:23 7:26 7:29 7:35 7:38 7:42 7:46 7:50 7:56 8:00

– 5 4 7 4 8 3 4 7 4 4 3 3 3 5 3 4 3 3 3 3 3 3 6 3 4 4 4 6 4

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

6:00 6:05 6:10 6:14 6:18 6:22 6:26 6:30 6:34 6:38 6:42 6:46 6:50 6:54 6:58 7:01 7:04 7:07 7:10 7:13 7:16 7:19 7:23 7:27 7:31 7:35 7:39 7:43 7:47 7:51 7:55 7:59

– 5 5 4 4 4 4 4 4 4 4 4 4 4 4 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

6:00 6:04 6:08 6:12 6:16 6:20 6:24 6:28 6:32 6:36 6:40 6:44 6:48 6:52 6:56 7:00 7:04 7:08 7:12 7:16 7:20 7:24 7:28 7:32 7:36 7:40 7:44 7:48 7:56 8:00

– 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 8 4

Table 5 Distribution of the 44,025 HAR Data. Workday

Each Inbound Bus Stop of Bus Route B1

Each Inbound Bus Stop of Bus Route B11

Each Inbound Bus Stop of Bus Route B23

Monday Tuesday Wednesday Thursday Friday Sum for one bus stop of one bus route Sum for all inbound stops of one bus route Total for all inbound stops of three bus routes

116 116 116 87 87 522

62 62 93 93 93 403

116 116 116 87 116 551

522 × 30 = 15660

403 × 28 = 11284

551 × 31 = 17081

x (5) (F ) = ξ0 +

44025

γ α [1−(1−F ) β]− [1−(1−F )−δ ] + η [1−(1−F )3] ,∀ F ∈ [0,1) δ β

where α,β,γ ,δ ,η are five parameters to be estimated using HAR data. Since function and −δ , without loss of generality, we assume that β ⩾ −δ , namely:

β+δ⩾0

(22)

x (5)

(F ) is unchanged by transformation between β (23)

As parameters β and δ are the dominators, we should impose that

β ≠ 0 and δ ≠ 0

(24)

Note that the conditions for parameters β and δ , as shown in Eqs. (23) and (24), are the same as the Wakeby distribution. To be a 230

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

Fig. 4. Histograms of the HAR Data at Each Inbound Bus Stop of Route B1 On Monday.

Table 6 Exemplify of the Fitting Results. No.

Number of estimated parameters

Distribution

Number of rejections

No.

Number of estimated parameters

Distribution

Number of rejections

1 2 3 4 5 6 7 8 9 10 11

2

Truncated Normal Log-Logistic Erlang Inverse Gaussian Fatigue Life Lognormal Pearson 5 Weibull Frechet Pert Gamma

28 30 29 25 30 30 30 21 30 17 13

12 13 14 15 16 17 18 19

3

Burr Johnson SB Generalized Gamma Beta Pearson 6 Kumaraswamy Phased Bi-Exponential Wakeby

10 11 19 15 22 14 28 6

4

Table 7 Number of Rejected GoF Tests for Wakeby Distribution With ξ = −1.

Monday Tuesday Wednesday Thursday Friday Average Standard Deviation

Route B1 (30 Stops)

Route B11 (28 Stops)

Route B23 (31 Stops)

6 (6/30 = 20.00%) 5(5/30 = 16.66%) 14 (14/30 = 46.67%) 16 (16/30 = 53.33%) 1 (1/30 = 3.33%) 8.4 (28.00%) 6.35 (31.16%)

3 (3/28 = 10.71%) 3 (3/28 = 10.71%) 3 (3/28 = 10.71%) 4 (4/28 = 14.29%) 3 (3/28 = 10.71%) 3.40 (12.14%) 0.55 (1.96%)

5 (5/31 = 16.67%) 7 (7/31 = 16.67%) 7 (7/31 = 16.67%) 2 (2/31 = 3.23%) 6 (6/31 = 19.35%) 5.40 (17.42%) 2.07 (6.69%)

valid quantile function, x (5) (F ) should be non-decreased with respect to F ∈ (0,1) . Nevertheless, the non-decreasing mononicity of function x (5) (F ) depends on the values taken by five parameters α,β,γ ,δ ,η . The following proposition gives a sufficient condition for these five parameters to guarantee the quantile functionality of x (5) (F ) with finite mean.

231

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

Table 8 Relations between PWMs and the Parameters of Tailored Wakeby-type Distributions: x (4) (F ) , x (3) (F ) and x (2) (F ) . Tailored Wakeby-type Distribution

β and −δ are the larger and smaller roots, respectively, for the quadratic equation:

x (4) (F )

(N2(4) C3(4)−N3(4) C2(4) ) z 2 + (N1(4) C3(4)−N3(4) C1(4) ) z + (N1(4) C2(4)−N2(4) C1(4) ) = 0 where

(72)

N1(4) = α 0(4)−16α1(4) + 27α2(4)−2ξ0 , N2(4) = α 0(4)−8α1(4) + 9α2(4) , N3(4) = α 0(4)−4α1(4) + 3α2(4) , C1(4) = 8α1(4)−54α2(4) + 64α3(4)−2ξ0 ,

C2(4) = 4α1(4)−18α2(4) + 16α3(4) , C3(4) = 2α1(4)−6α2(4) + 4α3(4)

Tailored Wakeby-type Distribution

x (3) (F )

γ = [(α 0(4)−ξ0)(1 + β )−(2α1(4)−ξ0)(2 + β )](1−δ )(2−δ )/(β + δ )

(73)

α = [(α 0(4)−ξ0)(1−δ )−(2α1(4)−ξ0)(2−δ )](1 + β )(2 + β )/(β + δ )

(74)

β = −(−2ξ0 + 4α 0(3)−40α1(3) + 54α2(3) )/(4α 0(3)−20α1(3) + 18α2(3) )

(75)

α = (ξ0 + 10α1(3)−18α2(3) )(2 + β )(3 + β )/(3−β ) η=

Tailored Wakeby-type Distribution

β=

x (2) (F )

α=

4 (−ξ0 3

+ α 0(3)−

α ) 1+β

(76)

(77)

−(α 0(2)−4α1(2) + ξ0)/(α 0(2)−2α1(2) ) (α 0(2)−ξ0)(1 + β ) (79)

(78)

Proposition 2. x (5) (F ) defined in Eq. (22) is the quantile function of a probability distribution with finite mean if

β⩾3

(25)

δ < 1 and δ ≠ 0

(26)

β+δ⩾0

(27)

α + γ + 3η ⩾ 0

(28)

γ⩾0

(29)

η⩾0

(30)

x (5)

(F ) mathematically meaningful because reciprocals of both Proof. Eqs. (25) and (26) imply that β ≠ 0 and δ ≠ 0 , which make parameters β and δ are involved in x (5) (F ) . The continuously differentiability of x (5) (F ) leads to the following sufficient and necessary conditions for x (5) (F ) to be the quantile function of a probability distribution: x (5) ′ (F ) =

d [x (5) (F )] ⩾ 0 ,∀ F ∈ (0,1) dF

(31)

It can be easily seen that

x (5) ′ (F )= α (1−F ) β − 1 + γ (1−F )−δ − 1 + 3η (1−F )2 = (1−F )−δ − 1 {(α + γ + 3η)(1−F ) β + δ + γ [1−(1−F ) β + δ ] + 3η [(1−F )δ + 3−(1−F ) β + δ ]}

(32)

Because (1−F )−δ − 1 > 0 for all F satisfying 0 < F < 1, Eq. (31) is thus equivalent to

(α + γ + 3η)(1−F ) β + δ + γ [1−(1−F ) β + δ ] + 3η [(1−F ) δ + 3−(1−F ) β + δ ] ⩾ 0 ,∀ F ∈ (0,1)

(33)

According to Eq. (27), we have

0 ⩽ (1−F ) β + δ ⩽ 1 ,∀ F ∈ (0,1)

(34)

0 ⩽ 1−(1−F ) β + δ ⩽ 1 ,∀ F ∈ (0,1)

(35)

In addition, for any β ⩾ 3, we also have

(1−F )δ + 3−(1−F ) β + δ ⩾ 0 ,∀ F ∈ (0,1) Eqs. (28)–(30) can, therefore, assure that x function x (5) (F ) can be calculated by 1

1

E (x (5) )= ∫0 x (5) (F ) dF = ∫0 {ξ + α

γ

α

(36) (5) ′

(F ) ⩾ 0 for any F ∈ (0,1) .The mean of random variable x (5) with the quantile

γ α [1−(1−F ) β]− δ [1−(1−F )−δ ] β γ

+ η [1−(1−F )3]} dF

η

ξ + β − δ + η + β ln(1−F )|01 − δ ln(1−F )|01 + 4 , β = −1 and δ = 1 ⎧ ⎪ γ η α γ α β = −1 and δ ≠ 1 ⎪ ξ + β − δ + η + β ln(1−F )|01 − δ (1 − δ ) (1−F )−δ + 1|01 + 4 , = γ η α γ α ⎨ ξ + β − δ + η + β (β + 1) (1−F ) β + 1|01 − δ ln(1−F )|01 + 4 , β ≠ −1 and δ = 1 ⎪ γ γ η α α 1 1 β δ + − + ⎪ξ + − + η + (1−F ) |01 − δ (1 − δ ) (1−F ) |01 + 4 , β ≠ −1 and δ ≠ 1 β δ β (β + 1) ⎩

(37)

Hence, Eqs. (25) and (26) also ensure the existence of E (x (5) ) , namely, the tailored Wakeby-type distribution has a finite 232

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

mean. □ The tailored Wakeby-type distribution with five parameters, as shown in Eq. (22), can be degraded into the following three tailored Wakeby-type distributions with four parameter, three parameters and two parameters by setting η = 0 , or γ = 0 , or γ = 0 and η = 0 , respectively. [Tailored Wakeby-type distribution with four parameters]:

x (4) (F ) = ξ0 +

γ α [1−(1−F ) β]− [1−(1−F )−δ ] ,∀ F ∈ [0,1) δ β

(38)

[Tailored Wakeby-type distribution with three parameters]:

x (3) (F ) = ξ0 +

α [1−(1−F ) β] + η [1−(1−F )3] ,∀ F ∈ [0,1) β

(39)

[Tailored Wakeby-type distribution with two parameters]:

x (2) (F ) = ξ0 +

α [1−(1−F ) β] ,∀ F ∈ [0,1),α ⩾ 0,β ≠ 0 β

(40)

x (4)

(F ) is identical to Wakeby distribution defined in Eq. (12). In other words, Eqs. (17)–(21) form a sufficient condition In reality of four parameters α,β,γ ,δ to ensure that x (4) (F ) is a valid quantile function of a probability distribution with limited mean. As for x (3) (F ) , we can derive the following proposition by means of the similar techniques used to prove Proposition 2. Proposition 3. x (3) (F ) defined in Eq. (39) is the quantile function of a probability distribution with finite mean if

β ≠ 0

(41)

β > −1

(42)

3η + α ⩾ 0

(43) (44)

α⩾0 It is noted that

x (2)

(F ) is the Pareto distribution. Moreover, the following conditions can sufficiently ensue that

x (2)

(F ) is a valid

Table 9 Estimated Parameter Estimates for the Probability Distributions of the Random HAR at Each Bus Stop Lying Along Bus Route B1 on Monday. Stop No.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

Wakeby Distribution

Tailored Wakeby-type Distribution

α

β

γ

δ

α

β

γ

δ

η

30.677 29.968 15.789 15.696 13.504 9.803 7.263 6.639 4.435 2.749 2.145 1.875 1.764 1.716 1.672 1.442 1.408 1.348 1.345 1.140 1.241 1.205 1.273 1.220 1.205 1.251 1.250 1.232 1.243 1.472

34.997 33.519 18.682 19.035 17.595 13.146 9.736 8.094 4.892 1.648 0.897 0.657 0.569 0.520 0.478 0.283 0.256 0.210 0.228 0.029 0.108 0.076 0.122 0.085 0.072 0.102 0.101 0.086 0.098 0.287

0.175 0.152 0.241 0.269 0.328 0.383 0.388 0.342 0.303 0.002 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0.257 0.317 0.174 0.139 0.103 0.037 0.067 0.120 0.156 0.980 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

32.686 30.975 25.603 30.504 21.581 13.077 8.116 6.193 4.435 2.749 0.893 1.435 1.417 1.056 1.021 1.146 1.276 1.379 1.407 2.088 1.711 2.248 1.926 1.838 1.709 1.709 1.578 1.975 1.843 1.283

38.453 35.254 43.327 55.695 41.551 28.608 20.529 13.083 4.892 1.648 0.288 0.472 0.420 0.205 0.163 0.143 0.196 0.224 0.256 0.397 0.304 0.465 0.379 0.334 0.281 0.291 0.240 0.376 0.339 0.202

0.158 0.142 0.099 0.102 0.167 0.201 0.218 0.220 0.303 0.002 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0.293 0.337 0.460 0.455 0.339 0.269 0.270 0.271 0.156 0.980 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0.049 0.028 0.444 0.493 0.444 0.498 0.556 0.503 0 0 0.584 0.208 0.169 0.337 0.338 0.161 0.072 −0.017 −0.033 −0.515 −0.257 −0.552 −0.349 −0.337 −0.279 −0.251 −0.182 −0.401 −0.326 0.102

Number of Rejected K-S Tests (SL = 0.05): 6

Number of Rejected K-S Tests (SL = 0.05): 4

233

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

quantile function of a probability distribution with finite mean:

As

β ≠ 0

(45)

α⩾0

(46)

β > −1

(47)

x (i)

(F ) (i = 2,3,4,5) are all valid quantile functions under the mild above-mentioned sufficient conditions, we have

x (i)

(F ) ⩾ ξ0 ,∀ F ∈ [0,1),i = 2,3,4,5

(48)

Remark 1. Tailored Wakeby-type Distributions are robust because of different parameters estimated. With various data, it is flexible to determine the number of parameters for the bus headway adherence issue for a given bus route with a number of bus stops. Moreover, according to the empirical experiments, more parameters adopted generally returns better fitting results (see experimental results that will be demonstrated in Tables 9–11). 4.1. Parameter estimation of the tailored Wakeby-type distributions The tailored Wakeby-type distributions defined above do not possess an explicit expression for their CDFs and PDFs except their analytical quantile functions. The L-moments based method, therefore, can be used to estimate those parameters involved in the tailored Wakeby-type distributions from sample data accordance to what Hosking (1990) has done for Wakeby distribution. Let us briefly introduce the L-moments of a random variable and its relationship with the probability weighted moments (Greenwood et al., 1979), and details can be found in Hosking and Wallis (1997). Let X be a real-valued random variable with CDF F (x ) and quantile function x (F ) , and X1: n ⩽ X2: n ⩽ …⩽Xn : n be the order statistics of a random sample of size n drawn from the distribution of random variable X . The L-moments of random variable X are defined by the quantities: r−1

λr = r −1 ∑ (−1) k k=0

(r−k 1) E (X

r − k : r ),r

= 1,2,… (49)

where Table 10 Estimated Parameter Estimates for the Probability Distributions of the Random HAR at Each Bus Stop Lying Along Bus Route B11 on Wednesday. Stop No.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

Wakeby Distribution

Tailored Wakeby-type Distribution

α

β

γ

δ

α

β

γ

δ

η

26.889 26.041 18.113 16.577 18.511 11.684 7.152 4.410 4.766 4.559 5.397 3.860 3.356 3.386 3.290 2.330 1.601 1.653 1.470 1.525 1.506 1.489 1.431 1.601 1.419 1.449 1.453 1.632

29.412 28.738 22.993 21.119 25.980 16.519 10.555 7.419 8.475 8.047 8.375 5.520 4.183 4.505 4.469 1.622 0.493 0.539 0.373 0.431 0.396 0.385 0.323 0.470 0.272 0.301 0.308 0.462

0.107 0.127 0.346 0.351 0.444 0.421 0.492 0.697 0.745 0.728 0.529 0.505 0.408 0.458 0.472 0.093 0 0 0 0 0 0 0 0 0 0 0 0

0.270 0.193 −0.213 −0.202 −0.221 −0.103 −0.134 −0.315 −0.348 −0.314 −0.083 −0.076 0.017 −0.007 −0.020 0.483 0 0 0 0 0 0 0 0 0 0 0 0

26.889 26.041 34.478 25.841 22.731 12.067 7.073 4.410 4.766 4.559 5.397 3.860 3.356 3.386 3.290 2.330 3.389 2.700 3.003 2.326 1.899 1.609 1.620 2.094 3.344 3.141 2.947 1.611

29.412 28.738 61.395 44.986 36.876 18.443 11.867 7.419 8.475 8.047 8.375 5.520 4.183 4.505 4.469 1.622 1.035 0.891 0.874 0.725 0.553 0.437 0.403 0.656 0.863 0.835 0.795 0.453

0.107 0.127 0.104 0.130 0.314 0.372 0.436 0.697 0.745 0.728 0.529 0.505 0.408 0.458 0.472 0.093 0 0 0 0 0 0 0 0 0 0 0 0

0.270 0.193 0.261 0.197 −0.072 −0.053 −0.085 −0.315 −0.348 −0.314 −0.083 −0.076 0.017 −0.007 −0.020 0.483 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0.462 0.424 0.208 0.099 0.135 0 0 0 0 0 0 0 0 0 −0.791 −0.471 −0.709 −0.377 −0.192 −0.060 −0.096 −0.234 −0.905 −0.797 −0.708 0.010

Number of Rejected K-S Tests (SL = 0.05): 11

Number of Rejected K-S Tests (SL = 0.05): 2

234

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

Table 11 Estimated Parameter Estimates for the Probability Distributions of the Random HAR at Each Bus Stop Lying Along Bus Route B23 on Friday. Stop No.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

Wakeby Distribution

Tailored Wakeby-type Distribution

α

β

γ

δ

α

β

γ

δ

η

24.368 58.358 29.803 13.873 8.642 7.990 7.904 7.379 7.268 6.908 6.594 5.268 5.113 3.209 2.837 2.682 1.917 1.719 1.737 1.717 1.753 1.814 1.588 1.699 1.694 1.736 1.705 1.752 1.777 1.861 2.004

25.742 72.695 37.070 17.198 10.898 11.413 11.209 11.581 11.982 11.621 11.542 8.536 7.940 4.008 3.181 4.279 2.258 0.626 0.618 0.596 0.627 1.351 2.958 0.573 0.558 0.626 0.596 0.639 0.662 0.741 0.878

0.054 0.334 0.255 0.297 0.315 0.438 0.436 0.564 0.631 0.658 0.705 0.696 0.629 0.412 0.352 0.657 0.572 0 0 0 0 0.254 0.785 0 0 0 0 0 0 0 0

0.400 −0.612 −0.175 −0.232 −0.114 −0.197 −0.200 −0.324 −0.393 −0.414 −0.441 −0.472 −0.378 −0.052 0.004 −0.244 −0.239 0 0 0 0 0.185 −0.150 0 0 0 0 0 0 0 0

26.941 195.362 29.803 17.429 8.642 7.990 7.904 7.543 7.950 7.525 7.615 6.159 5.356 3.209 2.837 2.682 1.917 3.922 1.895 1.759 1.446 1.814 1.588 1.564 1.594 1.089 1.123 0.973 0.922 0.689 0.475

30.508 279.046 37.070 28.867 10.898 11.413 11.209 15.906 20.122 20.089 23.310 28.116 23.802 4.008 3.181 4.279 2.258 1.237 0.680 0.613 0.494 1.351 2.958 0.516 0.517 0.312 0.316 0.247 0.225 0.089 −0.053

0.029 0.092 0.255 0.102 0.315 0.438 0.436 0.398 0.383 0.396 0.398 0.288 0.280 0.412 0.352 0.657 0.572 0 0 0 0 0.254 0.785 0 0 0 0 0 0 0 0

0.572 −0.008 −0.175 0.195 −0.114 −0.197 −0.200 −0.173 −0.172 −0.189 −0.186 −0.086 −0.036 −0.052 0.004 −0.244 −0.239 0 0 0 0 0.185 −0.150 0 0 0 0 0 0 0 0

0.106 0.280 0 0.391 0 0 0 0.303 0.413 0.430 0.488 0.732 0.723 0 0 0 0 −0.929 −0.072 −0.020 0.146 0 0 0.064 0.048 0.317 0.286 0.385 0.421 0.582 0.754

Number of Rejected K-S Tests (SL = 0.05): 6

E (Xr − k : r )=

r! (r − k − 1) ! k !

∫ {x [F (x )]r − k − 1 [1−F (x )]k } dF (x )

=

r! (r − k − 1) ! k !

1 ∫0 {x (F ) F r − k − 1 (1−F )k} dF

Number of Rejected K-S Tests (SL = 0.05): 3

(50)

The L-moments to describe and characterize a probability distribution is well justified by the following fundamental proposition (Hosking, 1990): Proposition 4. The L-moments λr (r = 1,2,…) exist if and only if random variable X has finite mean and they can fully characterize random variable X provided that random variable X has a finite mean. The L in “L-moments” emphasizes that λr is a linear function of the expected order statistics. The probability weighted movements (PWMs) of random variable X is defined by

αr = E [X {1−F (X )}r ] r = 0,1,…

(51)

The PWMs and L-moments have the linear relations; for example:

λ1 = α 0

(52)

λ2 = α 0−2α1

(53)

λ3 = α 0−6α1 + 6α2

(54)

λ 4 = α 0−12α1 + 30α2−20α3

(55)

λ5 = α 0−20α1 + 90α2−140α3 + 70α4

(56)

4.1.1. Relations between PWMs and parameters of the tailored Wakeby-type distributions We first derive the first five L-moments shown below for the tailored Wakeby-type distribution x (5) (F ) according to Eqs. (49), (50) and (52)–(56). 235

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al. 1

λ1(5)= E (X1:1) = ∫0 {x (5) (F )} dF = α 0(5) = ξ0 +

α (1 + β )

+

1

λ 2(5)= 2−1 ∑ (−1) k =

k=0 α (1 + β )(2 + β )

+

2

λ3(5)= 3−1 =

∑

(−1) k

γ (1 − δ )

3

+ 4η

(57)

(k1 ) E (X

2 − k :2 )

γ (1 − δ )(2 − δ )

+

(k2 ) E (X

k=0 α (1 − β ) (1 + β )(2 + β )(3 + β )

3 η 20

3 − k :3 )

+

= αo(5)−2α1(5) (58)

= α 0(5)−6α1(5) + 6α 2(5)

γ (1 + δ ) 1 − η (1 − δ )(2 − δ )(3 − δ ) 20

(59)

3

3 λ4(5)= 4−1 ∑ (−1) k ⎛ ⎞ E (X 4 − k :4 ) = α 0(5)−12α1(5) + 30α 2(5)−20α3(5) ⎝k ⎠ k=0 =

α (1 − β )(2 − β ) (1 + β )(2 + β )(3 + β )(4 + β ) 4

λ5(5) = 5−1

∑

(−1) k

k=0

=

+

γ (1 + δ )(2 + δ ) (1 − δ )(2 − δ )(3 − δ )(4 − δ )

(k4) E (X

5 − k :5 )

α (1 − β )(2 − β )(3 − β ) (1 + β )(2 + β )(3 + β )(4 + β )(5 + β )

+

+

1 η 140

(60)

= α 0(5)−20α1(5) + 90α 2(5)−140α3(5) + 70α4(5) γ (1 + δ )(2 + δ )(3 + δ ) (1 − δ )(2 − δ )(3 − δ )(4 − δ )(5 − δ )

(61)

α 0(5),α1(5),α 2(5),α3(5),α4(5)

where are the PWMs. Based on Eqs. (57)–(61), we can express parameters α,β,γ ,δ ,η as the functions of PWMs: α 0(5),α1(5),α 2(5),α3(5),α4(5) . Following the techniques of Hosking and Wallis (1997), we have the following proposition with the rigorous proof presented in Appendix A. Proposition 5. β and −δ are the larger and smaller roots, respectively, for the quadratic equation with one unknown z :

(N2(5) C3(5)−N3(5) C2(5) ) z 2 + (N1(5) C3(5)−N3(5) C1(5) ) z + (N1(5) C2(5)−N2(5) C1(5) ) = 0

(62)

N1(5) = 4α 0(5)−120α1(5) + 486α 2(5)−448α3(5) + 6ξ0

(63)

N2(5) = 4α 0(5)−60α1(5) + 162α 2(5)−112α3(5)

(64)

N3(5) = 4α 0(5)−30α1(5) + 54α 2(5)−28α3(5)

(65)

C1(5) = 40α1(5)−486α 2(5) + 1344α3(5)−1000α4(5) + 16ξ0

(66)

C2(5) = 20α1(5)−162α 2(5) + 336α3(5)−200α4(5)

(67)

C3(5) = 10α1(5)−54α 2(5) + 84α3(5)−40α4(5)

(68)

where

Having β and δ by solving the quadratic equation expressed by Eq. (62), we have

α=−

[(2ξ0−4α 0(5) + 40α1(5)−54α 2(5) )−(−4α 0(5) + 20α1(5)−18α 2(5) ) δ ](2 + β )(1 + β )(3 + β ) 2(β−3)(β + δ )

γ=

[(2ξ0−4α 0(5) + 40α1(5)−54α 2(5) ) + (−4α 0(5) + 20α1(5)−18α 2(5) ) β ](2−δ )(1−δ )(3−δ ) 2(−δ −3)(β + δ )

η=

4γ 1 4α (−4ξ0 + 4α 0(5)− − ) 3 1 + β 1−δ

(69)

(70)

(71)

(F ) , (F ) and We can also find the relations between PWMs and parameters of the tailored Wakeby-type distributions – x (2) (F ) – by the similar techniques used above for x (5) (F ) , which are presented in Table 8. Note that relations between PWMs and parameters of x (4) (F ) , i.e., Wakeby distribution, with ξ0 = 0 are derived by Hosking and Wallis (1997). x (4)

x (3)

4.1.2. Sample L-moments based parameter estimation method Let x1: n ⩽ x2: n ⩽ …⩽x n : n be an ordered random sample with size n drawn from the distribution of random variable X , the sample Lmoments denoted by lr (r ⩽ n ), are defined by 236

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

Fig. 5. L-moments Based Parameter Estimation Method.

237

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al. r−1

lr =

∑

(pr∗− 1,k × bk ),r = 1,2,…,n

(80)

k=0

where n

(i−1)(i−2)…(i−k ) bk = n−1 ∑ ⎡ x ⎤ ⎢ (n−1)(n−2)…(n−k ) i : n⎥ ⎦ i=1 ⎣

(81)

(kr ) ⎛⎝r +k k ⎞⎠

(82)

pr∗,k = (−1)r − k

Hosking (1990) demonstrates that the sample L-moments - lr ,r = 1,2,…,n - are the unbiased estimations of the L-moments: λr ,r = 1,2,…,n . Hence, we can employ the L-moments based method proposed by (Hosking, 1990) to estimate parameters of a proper tailored Wakeby-type distribution fitting the HAR data. The L-moments based parameter estimation method will first estimate five parameters - α,β,γ ,δ and η - for the tailored Wakebytype distribution with five parameters as follows: (i) calculate the first five sample L-moments for the HAR data according to Eq. (80); (ii) let each of the first five sample L-moments be a estimation of the corresponding L-moments, namely, λr̂ = lr ,r = 1,2,…,5; (iii) substitute λr (r = 1,2,…,5) in Eqs. (52)–(56) with λr̂ (r = 1,2,…,5) and solve the resultant equation to get the estimations of the 5 PWMs, ̂ (r = 1,2,…,5) (iv) substitute αr(5) (r = 1,2,…,5) in Eqs. (62)–(71) with αr(5) ̂ (r = 1,2,…,5) and obtain the estimations of the five αr(5) parameters denoted by α ,̂ β ,̂ γ ,̂ δ ̂ and η .̂ If these five estimated parameters satisfy the sufficient conditions in Eqs. (25)–(30), then stop. It is because the tailored Wakeby-type distribution with the estimated parameters is the probability distribution of the random HAR. Otherwise, we will perform the above procedure sequentially for the tailored Wakeby-type distribution with four parameters until to the two parameters. It is possible that none of all the Wakeby-type distributions can be a probability distribution of the random HAR. Details of the L-moments based parameter estimation method are illustrated in Fig. 5. 4.2. Numerical comparisons of the Wakeby distribution and the tailored Wakeby-type distribution To assess if the tailored Wakeby-type distributions with the parameters that are estimated by the L-moments based parameter estimation method can play a right role of the unified probability distribution for the random HAR, the three bus routes with 44,025 HAR data summarized in Table 5 will be used again. According to the four-step procedure shown in Fig. 2, we first estimate the parameters of the tailored Wakeby-type distributions with ξ0 = −1 by the L-moments based estimation method depicted in Fig. 5 using the 44,025 HAR data. Tables 9–11 respectively give the parameters estimated for bus route B1 on Monday, bus route B11 on Wednesday and bus route B23 on Friday for the illustration purpose. For comparing, the fitting results by the conventional Wakeby distributionare also shown in Tables 9–11. These tables indicate that the tailored Wakeby-type distributions can fit the HAR data remarkable better than the conventional Wakeby distribution in view of the number of rejected cases for the tailored Wakeby-type distributions. According to these results, the conventional Wakeby distributions with 2 or 4 parameters are replaced by the tailored Wakeby-type distributions with 3 or 5 parameters. K-S test results of all the cases (445) are shown in Table 12, which shows that the hypothesis of the Tailored Wakeby-type distribution is rejected for 42 cases by K-S statistics (SL = 0.05). Compared with 85 rejections for the conventional Wakeby distribution, 51% rejections have been declined. This significant improvement fully shows that the proposed Tailored Wakeby-type distribution can be redeemed as the unified probability distribution of the random HAR for each bus stop of a bus routes. 5. Applications of the random HAR with tailored Wakeby-type distributions The tailored Wakeby-type distribution can be regarded as an off-line probability distribution and captures the randomness of bus headways accurately. On one hand, using some tangible assessment indexes, it facilitates a transportation agency to enhance its performance of a bus service reliability. On the other hand, it could be connected to other transit modeling applications. That is, the Wakeby-type distribution can be incorporated into a transit simulation tool because it provides a probability description rather than the deterministic results of bus headways. Probability modeling of bus headways is more precise than deterministic models due to its uncertain nature. It thus has the potential to be used for the problems such as bus bunching, transit assignment and scheduling, particularly for bus online control. Following are examples about the application of the tailored wakeby-type distribution. Table 12 Number of Rejected K-S Tests for the Tailored Wakeby-type Distribution.

Monday Tuesday Wednesday Thursday Friday

Route B1 (30 Stops)

Route B11 (28 Stops)

Route B23 (31 Stops)

4 0 6 3 0

3 2 2 2 2

4 4 5 2 3

238

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

Given the tailored Wakeby-type distributions for the random HAR, probability of the HAR greater or lower than any given value x at a particular bus stop can be calculated by (83)

F (X ) = P (RHAR ⩽ x )

For a special given value x , two useful and novel indexes can be obtained by the derived probability distribution of the random HAR and they are elaborated below. Probability of Headway Adherence (PHA). The PHA is defined as probability that the HAR is kept within an acceptable range [x1,x2], where x1 and x2 are two pre-determined values. Ideally, bus headway is expected to be the same with scheduled headway. However, in a stochastic environment, deviation from schedule is unavoidable due to uncertain travel time and uneven passenger demand. It is thus reasonable to define an acceptable range within which HAR can be treated as reliable. The PHA can be calculated by (84)

PHA = F (x1)−F (x2) = P (x1 ⩽ RHAR ⩽ x2)

Probability of Bus Bunching (PBB). The PBB can be defined as the probability that two or more buses deployed by the same bus route arrive together or in close succession θ and it is useful in identifying those bus stops with frequent bunching phenomena. It is straightforward that parameter θ should be related to the scheduled headway. Hence, we define θ by (85)

θ = τH sch H sch

is the scheduled headway and τ ∈ (0,1) can be named as the coefficient of bus bunching. Value of τ conveys the tolerance where of bus bunching, and a smaller value of τ causes fewer trips to be sorted as bunching. Let H act denote actual bus headway, thus we have

PBB= P (0 ⩽ H act ⩽ θ) = P

(

−H sch H sch

⩽

H act − H sch H sch

⩽

τH sch − H sch H sch

)

= P (−1 ⩽ RHAR ⩽ τ −1) = F (τ −1)−F (−1) = F (τ −1)

(86)

We now analyze the bus bunching phenomena of bus route B1 introduced in Section 3.2 by using the PBB. After setting τ = 0.2, we have PBB = F(−0.8) according to Eq. (86). According to the tailored Wakeby-type distributions estimated for each bus stop of the bus route, we calculate PBB of each inbound bus stop on each work day, as shown in Fig. 6. According to Fig. 6, it can be clearly seen severity of the PBB at each bus stop by assuming the horizontal base line with PBB = 0.1 highlighted in the bold line. It provides a guideline of the bus stops at which bus bunching should be controlled. Fig. 6 also shows a general collar-shaped trend in terms of the bus bunching probability: at the beginning of B1 route, PBB is lower while it increases with every added bus stop; however, a slight fall around the 20th bus stop follows the sharp increase. The growth trend (the dashed line in Fig. 6) is coincident with the rule of thumb that bus bunching increases with the distance from the departure terminal since the variability may increase itself, which is also known as bunching effect: a lagging bus has to collect more passengers, and therefore tends to fall further behind (Daganzo and Pilachowski, 2011). While, an interesting finding is that the worsen does not continue until the last stop but recovers a bit at the end of the line, although there is no control point along the route. The recovery (the dotted line in Fig. 6) comes partly due to the general passenger pattern of bus route: less boarding at the end of route, which leads to less variance. Note that, B1 is a BRT route fully equipped with the dedicated bus lane, and the en-route travel reliability is thus mainly influenced by passenger demand. Another reason is that, as we mentioned before, derivers are expected to endeavor to maintain the scheduled arrival time at the end of each trip, which impels an improved schedule adherence. In addition to the application in bus bunching control point identification, PBB is also useful in bus bunching control strategies decision. Strategies for eliminating bus bunching can be classified into two groups: holding based on schedule and controlling based on real-time headway information. The former inserts slack time into schedule, and the slack is calculated so that buses can make up delay due to random recurrent disruptions (Zhao et al., 2006). If regular headway is the objective of control, early buses whose headway below a certain amount of time (e.g. the desired headway) are held for additional time at control points (Hickman, 2001).

Fig. 6. PBB at Each Inbound Bus Stop of Route B1 on Each Workday.

239

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

Although more slack time being planned leads to more reliable service, more severe in-vehicle delay and lower commercial speed would also be caused. To trade-off between service reliability and service efficiency, optimization on desired headway is the crux of the issue. There is no doubt that a real data derived probability distribution about headway, allowing analysts to see how often headways were shorter than any given threshold and how many trips would be held with a given desired headway, would provide meaning guideline on determination of the desired headway. The control strategy based on real-time headway information has been proposed by Xuan et al. (2011) and Daganzo (2009). Although adjusting is based on real-time information, the ideal achievable headway is not static and not even knowable in advance, it changes continually with traffic condition and passenger demand. It is thus reasonable to take the probability distribution of random HAR into consideration. 6. Conclusions and discussions A probabilistic distribution to describe the random HAR is necessary. As we stated in the Introduction (the second paragraph), bus headway is determined by a number of factors, and it is thus of great practical significance to investigate the bus headway adherence issue from a probabilistic perspective. However, in most cases, traditional distributions are unable to describe the bus headway in reality (see in the section of 1.1 Relevant Studies). Moreover, for the bus headway related studies, where need an assumption of headway distribution, a unified probabilistic distribution is necessary. This study addressed the bus headway adherence for bus routes from a probabilistic perspective owing to those uncertain parameters affecting bus headway adherence. After defining the random HAR, a customized four-step procedure to estimate the parameters of a probability distribution that can well describe the randomness of the bus HAR at each bus stop of the bus route was presented. According to our real case study for the three bus routes operated in Changzhou City of China with 44,025 HAR data in total, the 19 existing probability distributions were unable to meet our expectation. We thus proposed a tailored Wakeby-type distribution with five parameters by adding one more parameterized term in the conventional Wakeby distribution. The tailored Wakeby-type distribution with parameters estimated by the L-moments based parameter estimation method well fitted the HAR data, indicating that this five-parameter estimation is also applicable in practice with the L-moments based parameter estimation. Finally, applications of the tailored Wakeby-type distribution for qualifying the bus bunching phenomena were examined. In addition, a couple of aspects deserve to be highlighted: (i) To be exactly, a unified probabilistic distribution that covers different bus stops and bus route types is the result of this study according to the real data (rather than our hypothesis). Actually, in the initial stage of this research, as a rule of thumb, our assumption was that different route types may have different probabilistic distributions, even on different workdays, because of which we dealt with each type of route on each workday separately. But as a result of Section 3.2, this assumption is not true. And all trial results about these three routes suggested that the Wakeby distribution is relative superior to other traditional distributions in HAR fitting. As a result, we tailed such a Wakeby-type distribution. (ii) There is no guarantee that all HAR samples from different time periods of a bus stop could be described by the tailored Wakebytype distribution. Because there is much uncertainty in the parameters affecting bus headway adherence, and more uncertainty leads to more difficulty in distribution describing. Moreover, if data are available, then non-peak time data should also be tested to have a comparison with the peak time results. To be a matured probability distribution in HAR fitting, the proposed Wakeby-type distribution needs more works to be done: (i) Further research is expected on the transferability and generalizability of the proposed distribution. That is, case study in this paper should be continued for HAR data of different bus services, although the Wakeby-type distribution is hoped to be more powerful in distribution describing. (ii) Robustness of the proposed distribution is expected to research in depth. Method of parameter estimation has been provided in this paper, and factors affecting the values of the parameters and the level of GOF is hoped to be further studied. (iii) A substantive case study to show the application of Wakeby-type distribution can be conducted. Application of the probabilitybased measures in bus service diagnosing has shown the necessity of Wakeby-type distribution. Without such a distribution, these probability-based measures could not be applied. Because in most cases, traditional distributions are unable to describe the bus headway in reality, as we reviewed before. Therefore, a specific application of Wakeby-type distribution is desirable. Acknowledgements We would like to express our sincere gratitude to J.R.M. Hosking, the research staff member of IBM T.J. Watson Research Center, for proving us the valuable studies on the conventional Wakeby distribution analysis. We also thank the Changzhou Public Transportation Company for sharing us the AVL data of three bus routes. We are grateful to the three anonymous reviewers and the editor-in-chief for their valuables comments and suggestions made for the early versions of this study. The first author much appreciate the solid research supervision of the second author (Prof. Qiang Meng) when she conducted this study at National University of Singapore. This study is supported by the research project “Public Transit Bus and Driver Scheduling” (R-302-000-172-114).

240

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

Appendix A. Proof of Proposition 5 To express the five unknown parameters involved in Eqs. (57)–(61) by the first five L-moments, PWM defined in Eq. (51) could be employed. Substituting Eq. (22) into Eq. (51), we have 1

αr(5)= E [X (F ){1−F (X )}r ] = ∫0 X (F )[1−F ]r dF

(

)

1

α r+β+1

+

γ

α

α

1

= ξ0 + β − δ + η ∫0 (1−F )r dF − β ∫0 (1−F ) β + r dF + =

1 r+1

(ξ

0

+

γ r−δ+1

+

3η 3+r+1

γ δ

1 1 ∫0 (1−F )−δ + r dF −η ∫0 (1−F )3 + r dF

)

(A.1)

Then, we have

rαr(5) − 1 = ξ0 +

γ 3η α + + r+β r −δ 3+r

(A.2)

Multiply both sides of Eq. (A.2) by (r + β )(r −δ )(3 + r ) :

(rαr(5) − 1−ξ0 )(r + β )(r −δ )(3 + r ) = α (r −δ )(3 + r ) + γ (r + β )(3 + r ) + 3η (r + β )(r −δ )

(A.3)

Expand the brackets: 3 2 2 (rαr(5) − 1−ξ0 )[r + (β−δ + 3) r + (−βδ + 3β−3δ ) r −3βδ ] = (α + γ + 3η) r + [α (3−δ ) + γ (β + 3) + 3η (β−δ )] r −3αδ + 3γβ−3ηβδ

(A.4) And let:

m1 = β−δ + 3

(A.5)

m2 = −βδ + 3β−3δ

(A.6)

m3 = −3βδ

(A.7)

m4 = α + γ + 3η

(A.8)

m5 = α (3−δ ) + γ (β + 3) + 3η (β−δ )

(A.9)

m6 = −3αδ + 3γβ−3ηβδ

(A.10)

Substituting Eqs. (A.5)–(A.10) into Eq. (A.4), we have: (5) (5) (5) 2 2 r 3 (rαr(5) − 1−ξ0 ) + m1 r (rαr − 1−ξ ) + m2 r (rαr − 1−ξ ) + m3 (rαr − 1−ξ ) = m4 r + m5 r + m6

(A.11)

Write down the equation for r = 1,…,5:

(α 0(5)−ξ0) + m1 (α 0(5)−ξ0) + m2 (α 0(5)−ξ0) + m3 (α 0(5)−ξ0) = m4 + m5 + m6

(A.12)

8(2α1(5)−ξ0)

(A.13)

+

4m1 (2α1(5)−ξ0)

+

2m2 (2α1(5)−ξ0)

+

m3 (2α1(5)−ξ0)

= 4m4 + 2m5 + m6

27(3α 2(5)−ξ0) + 9m1 (3α 2(5)−ξ0) + 3m2 (3α 2(5)−ξ0) + m3 (3α 2(5)−ξ0) = 9m4 + 3m5 + m6

(A.14)

64(4α3(5)−ξ0) + 16m1 (4α3(5)−ξ0) + 4m2 (4α3(5)−ξ0) + m3 (4α3(5)−ξ0) = 16m4 + 4m5 + m6

(A.15)

125(5α4(5)−ξ0) + 25m1 (5α4(5)−ξ0) + 5m2 (5α4(5)−ξ0) + m3 (5α4(5)−ξ0) = 25m4 + 5m5 + m6

(A.16)

Subtract Eq. (A.13) from Eq. (A.12), Eq. (A.14) from Eq. (A.13), etc.:

α 0(5)−16α1(5) + 7ξ0 + m1 α 0(5)−8m1 α1(5) + 3m1 ξ0 + m2 α 0(5)−4m2 α1(5) + m2 ξ0 + m3 α 0(5)−2m3 α1(5) = −3m4−m5

(A.17)

16α1(5)−81α 2(5) + 19ξ0 + 8m1 α1(5)−27m1 α 2(5) + 5m1 ξ0 + 4m2 α1(5)−9m2 α 2(5) + m2 ξ0 + 2m3 α1(5)−3m3 α 2(5) = −5m4−m5

(A.18)

81α 2(5)−256α3(5)

(A.19)

+ 37ξ0 +

256α3(5)−625α4(5)

27m1 α 2(5)−64m1 α3(5)

+ 61ξ0 +

+ 7m1 ξ0 +

64m1 α3(5)−125m1 α4(5)

9m2 α 2(5)−16m2 α3(5)

+ 9m1 ξ0 +

+ m2 ξ0 +

16m2 α3(5)−25m2 α4(5)

3m3 α 2(5)−4m3 α3(5)

+ m2 ξ0 +

= −7m4−m5

4m3 α3(5)−5m3 α4(5)

= −9m4−m5

(A.20)

Subtract Eq. (A.18) from Eq. (A.17), Eq. (A.19) from (A.18), etc.:

α 0(5)−32α1(5) + 81α 2(5)−12ξ0 + m1 α 0(5)−16m1 α1(5) + 27m1 α 2(5)−2m1 ξ0 + m2 α 0(5)−8m2 α1(5) + 9m2 α 2(5) + m3 α 0(5)−4m3 α1(5) + 3m3 α 2(5) (A.21)

= 2m4

241

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

16α1(5)−162α 2(5) + 256α3(5)−18ξ0 + 8m1 α1(5)−54m1 α 2(5) + 64m1 α3(5)−2m1 ξ0 + 4m2 α1(5)−18m2 α 2(5) + 16m2 α3(5) + 2m3 α1(5)−6m3 α 2(5) + 4m3 α3(5) = 2m4

(A.22)

81α 2(5)−512α3(5) + 625α4(5)−34ξ0 + 27m1 α 2(5)−128m1 α3(5) + 125m1 α4(5)−2m1 ξ0 + 9m2 α 2(5)−32m2 α3(5) + 25m2 α4(5) + 3m3 α 2(5)−8m3 α3(5) + 5m3 α4(5) = 2m4

(A.23)

Subtract Eq. (A.22) from Eq. (A.21), Eq. (A.23) from Eq. (A.22):

α 0(5)−48α1(5) + 243α 2(5)−256α3(5) + 6ξ0 + m1 α 0(5)−24m1 α1(5) + 81m1 α 2(5)−64m1 α3(5) + m2 α 0(5)−12m2 α1(5) + 27m2 α 2(5)−16m2 α3(5) + m3 α 0(5)−6m3 α1(5) + 9m3 α 2(5)−4m3 α3(5) = 0

(A.24)

16α1(5)−243α 2(5) + 768α3(5)−625α4(5) + 16ξ0 + 8m1 α1(5)−81m1 α 2(5) + 192m1 α3(5)−125m1 α4(5) + 4m2 α1(5)−27m2 α 2(5) + 48m2 α3(5)−25m2 α4(5) + 2m3 α1(5)−9m3 α 2(5) + 12m3 α3(5)−5m3 α4(5) = 0

(A.25)

Substituting Eqs. (A.5)–(A.7) into Eqs. (A.24), (A.25) yields that:

(4α 0(5)−120α1(5) + 486α 2(5)−448α3(5) + 6ξ0) + (4α 0(5)−60α1(5) + 162α 2(5)−112α3(5) )(β−δ ) + (−2α 0(5)−30α1(5) + 54α 2(5)−28α3(5) )(−βδ ) = 0 (A.26)

(40α1(5)−486α 2(5) +

+

1344α3(5)−1000α4(5)

(10α1(5)−54α 2(5)

+

+ 19ξ0) +

84α3(5)−40α4(5) )(−βδ )

(20α1(5)−162α 2(5)

+

336α3(5)−200α4(5) )(β−δ )

=0

(A.27)

Eqs. (A.26), (A.27) can be rewritten as follows:

N1(5) + N2(5) (β−δ ) + N3(5) (−βδ ) = 0

(A.28)

C1(5) + C2(5) (β−δ ) + C3(5) (−βδ ) = 0

(A.29)

N1(5) = 6ξ0 + 4α 0(5)−120α1(5) + 486α 2(5)−448α3(5)

(A.30)

N2(5)

(A.31)

where

=

4α 0(5)−60α1(5)

+

162α 2(5)−112α3(5)

N3(5) = 4α 0(5)−30α1(5) + 54α 2(5)−28α3(5)

(A.32)

C1(5) = 16ξ0 + 40α1(5)−486α 2(5) + 1344α3(5)−1000α4(5)

(A.33)

C2(5) = 20α1(5)−162α 2(5) + 336α3(5)−200α4(5)

(A.34)

C3(5)

(A.35)

=

10α1(5)−54α 2(5)

+

84α3(5)−40α4(5)

Solve for β−δ and −βδ , obtaining

β−δ =

N3(5) C1(5)−N1(5) C3(5) N2(5) C3(5)−N3(5) C2(5)

β (−δ ) =

(A.36)

N1(5) C2(5)−N2(5) C1(5) N2(5) C3(5)−N3(5) C2(5)

(A.37)

Therefore, the quadratic equation

(N2(5) C3(5)−N3(5) C2(5) ) z 2 + (N1(5) C3(5)−N3(5) C1(5) ) z + (N1(5) C2(5)−N2(5) C1(5) ) = 0

(A.38)

has the right sides of Eqs. (A.36) and (A.37) as the sum and product of its roots. Thus, β and –δ are the roots of Eq. (A.38). Given β and δ, write down the expression in Eq. (A.2) for r = 1,2,3:

α 0(5) = ξ0 +

γ 3η α + + 1+β 1−δ 4

(A.39)

2α1(5) = ξ0 +

γ 3η α + + 2+β 2−δ 5

(A.40)

3α 2(5) = ξ0 +

γ 3η α + + 3+β 3−δ 6

(A.41)

Solve above equation set for α, γ and k, we have

242

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

α=−

γ=

η=

[(2ξ0−4α 0(5) + 40α1(5)−54α 2(5) )−(−4α 0(5) + 20α1(5)−18α 2(5) ) δ ](2 + β )(1 + β )(3 + β ) 2(β−3)(β + δ )

[(2ξ0−4α 0(5)

+

40α1(5)−54α 2(5) )

(−4α 0(5)

(A.42)

20α1(5)−18α 2(5) ) β ](2−δ )(1−δ )(3−δ )

+ + 2(−δ −3)(β + δ )

(A.43)

4γ ⎞ 1 ⎛ (5) 4α ⎜4α − ⎟ 0 −4ξ0− 3⎝ 1 + β 1−δ ⎠

(A.44)

References Adebisi, O., 1986. A mathematical model for headway variance of fixed-route buses. Transport. Res. Part B: Methodol. 20 (1), 59–70. Ap. Sorratini, J., Liu, R., Sinha, S., 2008. Assessing bus transport reliability using micro-simulation. Transport. Plan. Technol. 31 (3), 303–324. Barabino, B., Di Francesco, M., Mozzoni, S., 2013. Regularity diagnosis by automatic vehicle location raw data. Public Transp. 4 (3), 187–208. Bellei, G., Gkoumas, K., 2010. Transit vehicles' headway distribution and service irregularity. Public Transp. 2 (4), 269–289. Bian, B., Zhu, N., Ling, S., Ma, S., 2015. Bus service time estimation model for a curbside bus stop. Transport. Res. Part C: Emerg. Technol. 57, 103–121. Byon, Y.J., Cortés, C.E., Jeong, Y.S., Martínez, F.J., Munizaga, M.A., Zúñiga, M., 2017. Bunching and headway adherence approach to public transport with GPS. Int. J. Civil Eng. http://dx.doi.org/10.1007/s40999-017-0153-3. Ceder, A., 2007. Public Transit Planning and operation: Theory, Modelling and Practice. Butterworth-Heinemann, Burlington, Mass. Ceder, A., Golany, B., Tal, O., 2001. Creating bus timetables with maximal synchronization. Transp. Res. Part A 35 (10), 913–928. Chang, S.K.J., Hsu, S.C., 2004. Modeling of passenger waiting time in intermodal station with constrained capacity on intercity transit. Int. J. Urban Sci. 8 (1), 51–60. Changzhou Statistics, 2013, Changzhou National Economic and Social Development Statistics Bulletin 2012. < http://www.cztjj.gov.cn/node/tjsj_3/2013-3-1/ 13311515436710155.html > . Cortes, C.E., Burgos, V., Fernandez, R., 2010. Modelling passengers, buses and stops in traffic micro simulation: review and extensions. J. Adv. Transport. 44 (2), 72. Daganzo, C.F., 2009. A headway-based approach to eliminate bus bunching: systematic analysis and comparisons. Transport. Res. Part B: Methodol. 43 (10), 913–921. Daganzo, C.F., Pilachowski, J., 2011. Reducing bunching with bus-to-bus cooperation. Transport. Res. Part B: Methodol. 45 (1), 267–277. de Cea, J., Fernández, E., 1993. Transit assignment for congested public transport systems: an equilibrium model. Transport. Sci. 27 (2), 133–147. Dou, X.P., Meng, Q., Guo, X.C., 2015. Bus schedule coordination for the last train service in an intermodal bus-and-train transport network. Transp. Res. Part C 60, 360–376. Ehrlich, J.E., 2010. Applications of Automatic Vehicle Location Systems Towards Improving Service Reliability and Operations Planning in London, Massachusetts Institute of Technology. Estrada, M., Mensión, J., Aymamí, J.M., Torres, L., 2016. Bus control strategies in corridors with signalized intersections. Transp. Res. Part C 71, 500–520. Greenwood, J.A., Landwehr, J.M., Matalas, N.C., Wallis, J.R., 1979. Probability weighted moments: definition and relation to parameters of several distributions expressable in inverse form. Water Resour. Res. 15 (5), 1049–1054. Guo, G., Luo, H., Lin, X., and Feng, C., 2011. Headway-based evaluation of bus service reliability. In: 14th International IEEE Conference on Intelligent Transportation Systems, Washington, DC, USA, pp. 1864–1868. Hakkesteegt, P., Muller, Th.M.J., 1981. Research increasing regularity. Verkeerskundige werkdagen 415–436 (in Dutch). Hickman, M.D., 2001. An analytic stochastic model for the transit vehicle holding problem. Transport. Sci. 35 (3), 215–237. Hill, S.A., 2003. Numerical analysis of a time-headway bus route model. Phys. A 328 (1–2), 261–273. Hosking, J.R.M., 1986. The Theory of Probability Weighted Moments, Research Report RC 12210. IBM Research Division, Yorktown Heights, NY 10598. Hosking, J.R.M., 1990. L-Moments: analysis and estimation of distributions using linear combinations of order statistics. J. Roy. Stat. Soc.: Ser. B (Methodol.) 52 (1), 105–124. Hosking, J.R.M., Wallis, J.R., 1997. Regional Frequency Analysis: An Approach Based on L-Moments. Cambridge University Press, Cambridge. Houghton, J.C., 1978. Birth of a parent: the Wakeby distribution for modeling flood flows. Water Resour. Res. 14 (6), 1105. Ibarra-Rojas, O.J., Rios-Solis, Y.A., 2012. Synchronization of bus timetabling. Transport. Res. Part B: Methodol. 46 (5), 599–614. Kang, L., Meng, Q., 2017. Two-phase decomposition method for the last train departure time choice in subway networks. Transp. Res. Part B 104, 568–582. Land Transport Authority (LTA) of Singapore, 2013. Singapore Land Transport: Statistics in Brief 2013. < http://www.lta.gov.sg/content/dam/ltaweb/corp/ PublicationsResearch/files/FactsandFigures/Stats_in_Brief_2013.pdf > . Land Transport Authority (LTA) of Singapore, 2014. Annex A: presentation on BSRF Trial. < http://app.lta.gov.sg/data/apps/news/press/2014/20140124_BSRF (final2)-Annex.pdf > . Lin, J., Ruan, M., 2009. Probability-based bus headway regularity measure. Intell. Transp. Syst., IET 3 (4), 400–408. Liu, R., Sinha, S., 2007. Modelling urban bus service and passenger reliability. In: International Symposium on Transportation Network Reliability, The Hague, July 2007. McLeod, F., 2007. Estimating Bus passenger waiting times from incomplete bus arrivals data. J. Operat. Res. Soc. 58 (11), 1518–1525. Meng, Q., Qu, X., 2013. Bus dwell time estimation at bus bays: a probabilistic approach. Transport. Res. Part C: Emerg. Technol. 36, 61–71. Munoz, J.C., Cortes, C.E., Giesen, R., Sez, D., Delgado, F., Valencia, F., Cipriano, A., 2013. Comparison of dynamic control strategies for transit operations. Transport. Res. Part C: Emerg. Technol. 28, 101–113. Nagatani, T., 2001. Bunching transition in a time-headway model of a bus route. Phys. Rev. E 63 (3), 036115. Nagatani, T., 2003. Chaos and headway distribution of shuttle buses that pass each other freely. Phys. A 323, 686–694. Nesheli, M.M., Ceder, A., 2016. Use of real-time operational tactics to synchronize transfers in headway-based public transport service. Transp. Res. Rec. 2539. http:// dx.doi.org/10.3141/2539-12. Petersen, H.L., Larsen, A., Madsen, O., Petersen, B., Ropke, S., 2013. The simultaneous vehicle scheduling and passenger service problem. Transport. Sci. 47 (4), 603–616. Polus, A., 1975. Measurement of transportation system reliability: concepts and applications., PhD Dissertation, Northwestern University, Ann Arbor, pp. 252-252. Polus, A., 1979. An analysis of the headway distribution of an urban bus service. Traffic Eng. Control 20 (8/9), 419–421. Rashidi, S., Ranjitkar, P., 2013. Approximation and short-term prediction of bus dwell time. In: Conference of Eastern Asia Society for Transportation Studies, 9. Ross, S.M., 2009. Introduction to Probability and Statistics for Engineers and Scientists. Academic Press/Elsevier, Amsterdam. Saberi, M., Ali Zockaie, K., Feng, W., El-Geneidy, A., 2013. Definition and properties of alternative bus service reliability measures at the stop level. J. Public Transport. 16 (1), 97–122. Strathman, J.G., Dueker, K.J., Kimpel, T., 1999. Automated bus dispatching, operations control, and service reliability: baseline analysis. Transp. Res. Rec. 1666, 28–36. Szeto, W.Y., Jiang, Y., Wong, K.I., Solayappan, M., 2013. Reliability-based stochastic transit assignment with capacity constraints: formulation and solution method. Transp. Res. Part C 35, 286. Toledo, T., Cats, O., 2010. Mesoscopic simulation for transit operations. Transp. Res. Part C 18 (6), 896–908. Transportation Research Board, 2006. Transit Cooperative Research Program (TCRP) Report 113-Using Archived AVL-APC Data to Improve Transit Performance and Management. Transit Cooperative Research Program, Transportation Research Board, Washington, D.C. Transportation Research Board, 2013. Transit Cooperative Research Program (TCRP) Report 165-Transit Capacity and Quality of Service Manual (TCQSM), third ed.

243

Transportation Research Part C 86 (2018) 220–244

M. Zhang et al.

Transit Cooperative Research Program, Transportation Research Board, Washington, D.C. van Oort, N., 2011. Service reliability and urban public transport design, Delft University of Technology, the Netherlands. van Oort, N., van Nes, R., 2009. Regularity analysis for optimizing urban transit network design. Public Transp. 1 (2), 155–168. Washington, S., Karlaftis, M.G., Mannering, F.L., 2011. Statistical and Econometric Methods for Transportation Data Analysis. CRC Press, Boca Raton, Fla. Wu, J.H., Florian, M., Marcotte, P., 1994. Transit equilibrium assignment - a model and solution algorithms. Transport. Sci. 28 (3), 193–203. Wu, W., Liu, R., Jin, W., 2017. Modelling bus bunching and holding control with vehicle overtaking and distributed passenger boarding behaviour. Transp. Res. Part B 104, 175–197. Xuan, Y., Argote, J., Daganzo, C.F., 2011. Dynamic bus holding strategies for schedule reliability: optimal linear control and performance analysis. Transp. Res. Part B 45 (10), 1831–1845. Yan, S.Y., Tang, C.H., 2008. An integrated framework for intercity bus scheduling under stochastic bus travel times. Transport. Sci. 42 (3), 318–335. Yu, H., Chen, D., Wu, Z., Ma, X., and Wang, Y., 2016. Headway-based bus bunching prediction using transit smart card data. Transport. Res. Part C: Emerg. Technol., vol. 72, Supplement C, pp. 45–59. Zhao, J., Dessouky, M., Bukkapatnam, S., 2006. Optimal slack time for schedule-based transit operations. Transport. Sci. 40 (4), 529–539.

244