Origin-destination matrix estimation with a conditionally binomial model

We focus on a simple network abstraction to study how well the amount of local and through traffic can be inferred from minimum measurements of the total traffic. Our network, called “Minicity”, is shown in Fig. 1.

There are 3 sink/source nodes, denoted by W=West, C = Centre, E=East, and 4 measurement points, denoted by 1, 2, 3, 4. There are 6 flows with O-D pairs (W,C),(W,E),(C,E),(E,C),(E,W), and (C,W). At each measurement point, we observe the number of vehicles during fixed time intervals. In the West-East direction, denote the traffic flows F_(W,C),F_(C,E),F_(W,E) as X,Y,Z, and the observations as O₁=X+Z, O₂=Y+Z. The Minicity model is not observable (see [5] and [25]): the flows X,Y,Z and X−δ,Y−δ,Z+δ, with δ some constant, yield identical observations O₁ and O₂. However, the analysis is possible as a stochastic demand estimation problem, in which a traffic model has an important role in identifiability.

We focus on the flow structure of regular traffic consisting of a rather stable amount of vehicles. Assume that the observation data consist of a long sequence of measurements on usual working days. Then the traffic follows steady patterns at large, but stochastic variation makes each day a bit different. To get an idea of this variation, Fig. 2 presents the variation of two very illustrative quantities on working day Mondays at the resolution of 15 minutes. The upper plot shows the index of dispersion of counts (IDC; variance over mean) and the lower plot the coefficient of variation (CoV; standard deviation over mean).

First, recall that the IDC of a Poisson distribution always equals one. The upper plot shows that the traffic can be considered Poissonian during the most silent period of the night, but not at any other time. However, the IDC is relatively low (about 2) and stable also during the day between 9:30 and 15:30, whereas it is notably higher during the morning and afternoon rush hours. The non-Poissonian character of vehicular traffic is generally recognised, but the upper plot of Fig. 2 offers a more detailed picture. Although a popular extension of Poissonian modelling has been to introduce additional variables, ε and θ, to scale the relation between the variance and the mean, e.g., Var(X)=ε(E(X))^θ [7], we study a more specific model.

In a normal working day, the majority of the morning rush hour traffic consists of people’s regular drives to workplaces, educational institutions etc. There is a steady population of the same vehicles that create this traffic, but the actual presence of each particular vehicle can be considered as random. This suggests a binomial type of model; from n vehicles potentially appearing in a rush time interval I, each one actually appears with some probability p. Such a model is, however, immediately refuted by statistical data. The IDC of Bin(n,p) is 1−p, i.e., still lower than that of a Poisson distribution. Obvious shortcomings in the motivation of the above binomial model are that the probability of using a car (i) is not homogeneous in the population, and (ii) varies according to time by common external factors like weather, seasonal holiday activity, flu epidemics etc. The inhomogeneity challenge is serious, but it is not in the scope of this paper. Instead, we focus on the time variation by considering a doubly stochastic, conditionally binomial model, in which the parameter p is the same for all flows but varies from day to day as a stochastic process.

The presence of a common “activity factor” of a day is suggested by the data; we observed a strong positive correlation of measured vehicle counts at separate locations that cannot logically contain more than a negligible number of the same cars, i.e., at the opposite directions of the same road. For a comprehensive statistical analysis of the data utilised in this paper, see [16]. In particular, the correlation matrices of the measurement locations O₁,…,O₄ on Monday – Thursday mornings at 7:45 – 8:45 A.M. and 8:45 – 9:45 A.M.,

and

show the highest correlations between physically connected locations (motorway passing the city). However, all location pairs indicate considerable positive correlation. We interpret this common correlatedness to be due to a kind of “activity factor”, see also [2].

The assumption of an “activity factor” is also compatible with the lower plot of Fig. 2. The CoV is roughly constant throughout the working time, the rush hours included, despite of strongly differing traffic volumes. This means that the random variation of the daily sequence of traffic counts in a fixed quarter-hour is dominated by a common random factor, not by individual level variation, whose aggregated variance scales with n, not n². A detailed analysis of the variance in the conditionally binomial model is given in the next subsections.

Consider a fixed time interval I of a day, and the daily traffic counts in I at four measurement stations as depicted in Fig. 1. Assume that for each day k, k=1,2,…, there is a random variable γ_k with values in (0,1), interpreted as the “activity level” of day k. The sequence (γ_k) is assumed stationary and ergodic, i.e., its time average agrees with its average over the probability space. Let us then assume that, conditioned on γ_k, the daily West-East traffic flows in period I, denoted as \(X_{k}^{WE}\), \(Y_{k}^{WE}\), and \(Z_{k}^{WE}\) (see section 2.1) are independent and have distributions \(Bin(n_{X}^{WE},\gamma _{k})\), \(Bin(n_{Y}^{WE},\gamma _{k})\), \(Bin(n_{Z}^{WE},\gamma _{k})\), respectively, and similarly for the East-West counterparts. Assume Minicity to be uncongested so that the traffic via the city centre is not a realistic option for the through traffic and there is only one path for each O-D flow (see [11] for estimation under congestion).

To simplify, we use the notation \(n_{X} = n_{X}^{WE}, n_{Y} = n_{Y}^{WE}, n_{Z} = n_{Z}^{WE}, m_{X}=n_{X}^{EW}\), \(m_{Y}=n_{Y}^{EW}\), \(m_{Z}=n_{Z}^{EW}\). When considering the model of a single day, the index k is omitted. The n and m parameters represent the sizes of relatively stable populations of vehicles that can be used to make a trip from an origin to a destination, hence passing through some of the measurement points. The numbers n_X,n_Y,n_Z,m_X,m_Y,m_Z as well as the distribution of γ are considered as unknown.

One obvious candidate for modelling γ is the family of beta distributions, resulting in the so-called beta-binomial distributions. In a beta-binomial distribution with parameters (α,β,n) with large n, the squared CoV is approximately β/α, i.e. independent of n, which coincides well with the properties of the measured traffic.

Consider the measured quantities O₁=X+Z and O₂=Y+Z, and denote their first and second moments as m₁=E(O₁), m₂=E(O₂), v₁=Var(O₁), v₂=Var(O₂) and c₁₂=Cov(O₁,O₂). Using the equations

that hold for any square integrable U, V and any conditioning random variable \(\mathcal {A}\), we obtain the following expressions in the conditionally binomial model:

Corresponding equations hold for the East-West direction, we denote them by (2) ^EW. The role of covariance c₁₂ is essential in identifying n_Z as the other equations involve either n_X+n_Y or n_Y+n_Z. Note that because the value of γ is common for both directions, (2) and (2) ^EW together provide 10 equations for 8 unknowns: n_X,n_Y,n_Z,m_X,m_Y,m_Z,E(γ), and Var(γ).

Let us denote the IDCs and squared CoVs of O₁ and O₂, respectively, as

Now, straightforward computations yield the following result.

It is not hard to see that a similar model with conditionally Poissonian distributions (instead of binomial distributions) is not identifiable. We summarise briefly other model candidates and their identifiability in our context in Table 1.

We close this section by discussing the use of the conditionally binomial model in more complex networks. As an example, consider a line topology with 4 nodes {1,2,3,4}, 3 directed links {ℓ₁,ℓ₂,ℓ₃}:={(1,2),(2,3),(3,4)} and 6 routes \(\mathcal {R}:=\left \{{r_{1},r_{2},r_{3},r_{4},r_{5},r_{6},}\right \}:=\left \{ \{\ell _{1}\},\{\ell _{2}\},\{\ell _{3}\},\{\ell _{1},\ell _{2}\},\{\ell _{2},\ell _{3}\},\{\ell _{1},\ell _{2},\ell _{3}\}\right \}\). Line topologies are relevant for vehicular traffic as well as for train passenger traffic.

Assume that, conditionally on a random parameter γ, the amount of vehicles on route r_i in a time period is a random variable X_i with distribution Bin(n_i,γ), the X_is being independent given γ. Assume that the traffic O_i on link ℓ_i is measured, i=1,2,3. Denote by A the 3×6 matrix (the route incidence matrix) The counterpart of the moment Eqs. 2 can then be written in matrix form as

where O=(O₁,O₂,O₃)^T and n=(n₁,…,n₆)^T. Since Cov(O) is symmetric, (5) amounts to 9 equations for the 8 unknowns E(γ), Var(γ), n₁,…,n₆. Substituting the first equation into the last term of the second, (5) is transformed into

In fact, Eq. 6 is valid for any network and route system with conditionally binomial traffic. The (i,j)-element of A diag(n)A^T is

It is straightforward to check (for example, by computing the determinant) that the linear map (n₁,…,n₆)↦A diag(n)A^T is bijective. Thus, the vector n can be solved from the second equation of (6) as a linear combination of the elements of the matrix Feeding this expression of n into the first equation of (6) yields three (effectively, two) non-linear equations for the two first moments of γ. Thus, the model is identifiable (except for some singular cases, cf. Proposition 1). An analytical solution of the equations is probably hard to find, but the optimisation approach of section 3.2 can be applied in practical cases. We point out that the above solution works in two steps separating the vector n and γ. This topic will be elaborated in section 4.3.

Despite of the existence of the explicit solution (3), the model presents a principal challenge by depending on the accurate estimation of variances. To consider this in some detail, let γ₁,…,γ_N be independent copies of a random variable γ with values in (0,1), and for k=1,…,N, let U_k be an independent random variable with distribution Bin(n,γ_k). The unbiased sample variance has variance

where σ²=Var(U) and μ₄=E((U−EU)⁴) (e.g., [8]).

With N large and n not too small, the standard deviation of \(s^{2}_{N}\) is well approximated by

Although the value of p can thus be consistently estimated as \(p\approx 1-s^{2}_{N}/\bar {U}\), the simultaneous estimation of both parameters of the binomial distribution is known to be difficult. [10] showed that there is no unbiased estimator of either parameter alone. The difficulty is intuitively obvious for the case of small p and large n, because Bin(n,p) is then very close to Poisson(np), a distribution with a single parameter. However, the usually powerful principle of maximum likelihood fails in this problem also for larger values of p, because the likelihood function turns out to be almost constant on the set \(\left \{{(n,p)} : \: {np=\bar {U}}\right \}\).

When Var(γ)>0, an estimate of \(STD(s^{2}_{N})\) can be obtained by approximating U by a normal distribution with the same mean and variance. By Cochran’s theorem [9], N i.i.d. samples from the normal distribution satisfy

where \(\chi _{N-1}^{2}\) is the Chi-squared distribution with N−1 degreees of freedom. Figure 3 illustrates the confidence intervals of \(s^{2}_{N}/\sigma ^{2}\) at the risk levels α=0.05,0.1,0.2. Beyond N=10⁴, increasing the sample size decreases the uncertainty of the sample variance estimate extremely slowly.

The most frequently appearing element in the solution (3) is the difference ζ₁−ζ₂ of the relative variances of O₁ and O₂. In order to provide satisfactory inference, the errors of the estimates ζ₁ and ζ₂ should be at the level of a fraction of their difference, say, at most ξ|ζ₁−ζ₂|, with some not too high ξ∈(0,1).

Assume now that we have an i.i.d. N-sample from the conditionally binomial model, and let \(s_{N}^{2}\) be the unbiased sample variance of O₁. The requirement STD(ζ₁)≤ξ|ζ₁−ζ₂| is roughly equivalent to

Inserting to (10) σ²=(n_X+n_Z)²Var(γ)+(n_X+n_Z)(E(γ)−E(γ²)) and \(\text {Var}\left (\chi _{1}^{2}\right)=2\), writing with β_X+β_Y+β_Z=1 and simplifying, we obtain from (11) the condition for N:

This expression reveals an important feature of the conditionally binomial model. When

the required number N of samples does not depend heavily on n. In the opposite case, however, the required N grows quadratically in n. To illustrate magnitudes, assume that E(γ)=0.8 and STD(γ)=0.05 — then (E(γ)−E(γ²))/Var(γ)=63. Thus, the system size n affects the estimation precision adversely. Intuitively, the effect of the binomial fluctuations, on which the identification of n_Z is based, becomes with large n negligible in comparison to the effect of the variation of γ.

The estimation challenge is illustrated in Fig. 4. We simulated the model with γ being uniform on the interval (0.6,0.8) and the relative population sizes being (β_X,β_Y,β_Z)=(2/6,1/6,3/6), and studied the accuracy of estimation from (3) with system sizes n=60 and 600, which are below and above the critical magnitude given in (13). The model parameters were estimated in both cases in hundred samples of size N=5000, computing the estimates from increasing subsamples to see the speed of convergence. Figure 4 presents the mean relative errors of the estimates of E(γ), n=n_X+n_Y+n_Z and β_Z=n_Z/(n_X+n_Y+n_Z). The experiment suggests that the smaller system can be identified rather well with about 600 samples, whereas estimates from the larger system can easily be nonsense (e.g., negative) even when there are several thousands of samples (individual simulation runs were quite different, and the shapes of the point clouds varied). Note that the number of daily samples with some degree of homogeneity can hardly exceed 1000 in the real world.

We discuss the optimisation approach first in the context of the Minicity problem as it illustrates the main elements of the approach and serves our real traffic case study.

A natural approach to solve the traffic matrix estimation problem in both directions simultaneously is to minimise the squarred error of moment Eqs. 2 and (2) ^EW. However, the covariance equations pose an additional challenge in the minimisation approach: if the last equation in (2) produces an error ℓ=c₁₂−((n_X+n_Z)(n_Y+n_Z)Var(γ)+n_Z(E(γ)−E(γ²)) with parameters (n_X,n_Y,n_Z,E(γ),Var(γ)), then the modified parameters (n_X+ε,n_Y+ε,n_Z−ε,E(γ),Var(γ)), where ε=2ℓ/(E(γ)−E(γ²)), produce an equal error with an opposite sign. Thus the minimisation of the moment equations cannot uniquely determine optimal parameter values as the equations allow shifting traffic between local and through traffic.

To avoid this difficulty, we require that the covariance equations (the last equations of (2) and (2) ^EW) be solved exactly. In the general form, the problem is

where i,j run over traffic count measurement locations, \(\bar {O_{i}}\) is the sample mean over the days, and s² and q_ij stand for sample variance and covariance, respectively. Expressions containing the parameters are E(O_i),Var(O_i), and Cov(O_i,O_j), given by (2) and its counterpart (2) ^EW; κ is a weight and scaling parameter.

In our case, the solution of the quadratic equation for n_Z is

and a similar one holds for the East-West counterpart m_Z. The quadratic equation has exactly one positive solution, because c₁₂−Var(γ)n_Xn_Y>0 by the positivity of n_Z. To summarise, in the bidirectional traffic matrix estimation we search for (n_X,n_Y,m_X,m_Y,E(γ),Var(γ)) that minimise the cost function

where n_Z is given by (14) and m_Z respectively. As the estimation of the first moments is more robust than the estimation of the second moments, we require that the first moments to be fitted more accurately than the second moments. This is achieved by the weight parameter κ that multiplies the fitting error of the second order moments, see [12] and [13]. Additionally, κ scales the different orders of magnitude in the first and second moments.

The optimisation approach of a more general topology follows from Eq. 6 in section 2.3. The first task is to solve the vector n as explained in that section. The optimisation is applied to the observed link means to solve for E(γ) and Var(γ). Note that this process naturally separates the parameters expressing the number of potential vehicles in each OD pair and the common intensity variable γ. The discussion of our case study in section 4.3 elaborates more the interplay of these parameters, as well as benefits of separating them.

We studied the O-D estimation from a 15 minute interval traffic count data collected in the city of Tampere during 2011-2014. The city topology and the measurement points are illustrated in Fig. 5.

Unfortunately, various construction works caused interruptions in data and we selected the northern motorway with measurement locations Santalahti and Tampella for the traffic matrix estimation in the Minicity model. This road is a fast way to enter the city centre, but there is a large amount of traffic passing by the city centre.

An earlier study of the data suggested that the weekday morning traffic between 06:00 and 10:00 A.M. indicates the traffic flows most clearly. In the afternoons between 3:00 and 9:00 P.M. there can be some correlation between measurement points due to the traffic flows, but it is hard to detect from the general simultaneous activity in all the directions of the traffic. For other time periods, traffic bursts at the measurement locations near the city centre can be considered independent.

For the Minicity traffic matrix estimation, we select two morning periods 7:45 - 8:45 A.M. and 8:45 - 9:45 A.M. on weekdays from Monday through Thursday. The morning rush hour tends to shift slightly later on Fridays. Also, we exclude all public holidays, school holidays, isolated working days next to a public holiday, days between Christmas Eve and the New Year as well as months June – August due to lower traffic volumes. In this way, we end up with N=528 and N=522 days of traffic count measurements for the former and the latter rush hour period, respectively.

Although we removed atypical working days a priori there have been accidents and other events that produce outlier traffic counts that we removed by FAST-MCD covariance estimation [19]. We can safely assume that our data set contains less than 25% of contamination and configure the algorithm as recommended in [19]. To minimise the cost function of (15) with κ=0.00001, a global optimisation routine NMinimize in Mathematica is utilised in algorithm autoselection mode. Also, the region of minimisation is constrained to be (n_X,n_Y,m_X,m_Y,E(γ),Var(γ)) such that n_Z>0,n_Y>0,m_X>0,m_Y>0,E(γ)≤1,Var(γ)≥0,c₁₂−n_Xn_YVar(γ)>0,Var(γ)≥0,c₃₄−m_Xm_YVar(γ)>0. The two last constraints assure that n_Z and m_Z have positive solutions. The optimisation estimates the mean and the variance of the random variable γ. When plotting solutions against the measured data, we model γ by a beta distribution. The beta-binomial distribution has properties that fit well with our observations on the index of dispersion and the coefficient of the variation of the measured traffic, see section 2.2. However, some other distribution for γ may provide a better fit to the data.

We examine the robustness of the estimation results in two ways. First, we apply resampling to the original data, and then we examine the optimal parameter values with generated data.

Resampling: The resampling data are generated by selecting randomly without replacement 90% of the Monday – Thursday morning data to produce a collection of 1000 estimation samples. For morning traffic between 7:45 – 8:45 A.M., we study closely 162 estimation results with good fit to the data. In these results, the value of E(γ) varies from 0.55 to 0.95, i.e., the probability of executing a journey could also be clearly lower than the one estimated from the original data. Clustering analysis of the results indicates that most solutions are around the original data estimate, but there could be an alternative model candidate that has a lower E(γ)≈0.65. We note that this value of E(γ) would be closer to the results of the traffic matrix estimation for 8:45 – 9:45 A.M.. The traffic matrix probably evolves during the morning hours. For the peak traffic, the value of E(γ) is high, and when the traffic volumes decrease, also the value of E(γ) decreases.

Generated data: Next, we generate from a beta-binomial distribution 1000 samples each of size 500, with parameters that resemble the Tampere city case, and examine the accuracy of our traffic matrix estimation technique. Note that this realistic experiment is a challenging one; the parameters are not in the range of efficient inference identified in section 3.1, because the m_X and m_Y values for East–West local traffic are close to each other. Results of section 3.1 and parameters at the top line of Table 2 allow to calculate lower bounds for the required number of samples: in the West–East direction, Eq. 12 with ξ=1 gives rise to a lower bound of 100 119 samples. In comparison, in the East–West direction, the corresponding lower bound is 213 814 samples, which reflects the difficulty of the estimation when m_X and m_Y values are close to each other. Eq. 13 provides an upper bound for the total number of vehicles per traffic direction, n_X+n_Y+n_Z or m_X+m_Y+m_Z, so that the required number of samples does not grow quadratically in the total number of vehicles. These upper bounds are around 50 vehicles, whereas the estimated values for the first morning period give rise to 2811 vehicles in West-East direction and 1583 in the opposite direction. Thus this experiment is extremely challenging for the optimisation.

In the estimation experiment with simulated samples and known parameter values, the true value of the West–East ratio n_X/(n_X+n_Z) is captured very well, but the East–West ratio m_X/(m_X+m_Z) tends to be estimated too low. Further, the estimation tends to favour high n or m values together with low E(γ) values. In the model, the errors in the n and m parameters versus E(γ) compensate each other, which explains the good fits to the generated traffic observations. We conclude that the sample size of 500 is not large enough to estimate the exact parameters values in a robust and reliable manner.

In practice, there needs to be additional information to judge the correct levels of the n or m parameters and E(γ). There is freedom to increase one and to decrease the other while keeping the product intact, which is the basic challenge of Bin(n,p) estimation with both p and n unknown. [10] note that the estimation of n and p are linked together. On the other hand, in our model the variable γ, common to all flows, limits the degree of freedom. Our further experiments indicated that the balancing the overestimation/underestimation of the n or m values with E(γ) vanishes when some parameters are fixed to their true values.