Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
69,872 characters · 11 sections · 68 citation commands
Estimation of Recursive Route Choice Models with Incomplete Trip Observations
\newcolumntype{C}{>{\arraybackslash}p{4em}}
\affil[1]{\it School of Computing and Information Systems, Singapore Management University } \affil[2]{\it School of Computing, National University of Singapore } \affil[3]{\it Department of Civil & Mineral Engineering, University of Toronto }
{\bf Keywords:} Incomplete observation, recursive logit, nested recursive logit, expectation-maximization, decomposition-composition.
Travel demand modeling plays a vital role in transportation planning, operations, and management. In travel demand modeling, route choice is a process that associates the behaviors of selecting a route to travel from an origin to a destination. Such a route choice model can be used to assess travelers’ preferences regarding characteristics of different routes and travellers (e.g. travel time, travel cost, number of crossings, or age, gender, and trip purpose of the travellers). Route choice models can also be used to predict traffic flows and can have applications in traffic simulation osorio2015computationally or network pricing zimmermann2021strategic. The parameters of a route choice model are often estimated using path observations collected by, for instance, Global Positioning System (GPS) devices equipped on vehicles. It is commonly known that such trip observations would be incomplete due to several reasons, for example, the GPS signals are intermittent or some parts of the GPS data are censored from the analyst due to privacy concerns. This missing data issue is critical for the estimation of route choice models but to the best of our knowledge, this issue has not been properly investigated in the context of route choice modeling. We address the issue in this paper.
We focus on link-based recursive route choice models FosgFrejKarl13,MaiFosFre15 due to their several advantages, namely, the estimation does not require sampling of route choice alternatives and the models can be consistently estimated and easy to predict. Such a recursive model is based on the Rust's dynamic discrete choice framework RUST1987 where a path choice is modelled as a sequence of link choices and the parameters are estimated by combining maximum likelihood estimation (MLE) and dynamic programming zimmermann2020tutorial. A naive approach to deal with the missing-data issue under recursive models is to ignore the unconnected link observations and solve the MLE problem only based on connected ones. This approach can be done straightforwardly using existing estimation algorithms FosgFrejKarl13,MaiFosFre15,MaiBasFre15_DeC. However, unconnected segments in the observations would contain valuable information about traveller route choice behavior and would be important for the estimation. In this paper, we aim at developing new estimation algorithms that allow us to efficiently and practically recover information from these incomplete segments.
\noindentOur contributions: We make a number of contributions in this work. We first exploit the expectation-maximization (EM) method to deal with the missing-data issue. This method has been widely used to estimate statistical models with unobserved variables dempster1977maximum,mclachlan2019finite and we show how this scheme can be used in the context of route choice modeling under recursive models. More precisely, we design an EM algorithm that alternatively performs two steps where the first step is to compute an expected log-likelihood function by sampling paths to fill unconnected segments in the observations, and the second step maximizes that expected log-likelihood function. We further show that, under the recursive logit model FosgFrejKarl13, the expected log-likelihood function is concave, thus the second step can be done efficiently by convex optimization.
Moreover, observing that the EM algorithm would be expensive to perform as it requires solving several MLE problems, we further propose a new method to directly compute the probabilities of unconnected path segments. That is, we show that the probability of unconnected links can be obtained by solving one system of linear equations. This result allows us to compute the likelihood function of the incomplete observations where the missing segments are explicitly taken into consideration. Furthermore, to reduce the number of systems of linear equations to be solved and speed up the estimation process, we develop a new algorithm, called as the decomposition-composition (DC) method. The idea is to decompose each linear system above into two parts and compose the common parts to reduce the number of linear systems to be solved. The main advantage of our DC algorithm is that we only need to solve only one system of linear equations to obtain all the choice probabilities of unconnected segments on trips that share the same destination. We also show that such a system of linear questions always has a unique solution.
We test our proposed algorithms with two popular recursive route choice models, i.e., the recursive logit (RL) FosgFrejKarl13 and nested recursive logit (NRL) models, using a dataset collected from a real network. Our numerical experiments clearly show that our DC algorithm outperforms the EM and other baseline approaches in recovering missing information from observations. Our numerical results are based on the RL and NRL models, but our method can be applied to other existing recursive models in the literature, e.g., the recursive network multivariate extreme value Mai_RNMEV, discounted RL oyama2017discounted, or stochastic time-dependent RL mai2021RL_STD.
Literature reviews: The literature on route choice modeling covers path-based and link-based recursive models. Path-based models Prat09 rely on sampling of paths, thus the parameter estimates would depend on the sampling methods. Even when corrections are added to the choice probabilities to achieve consistent estimates, the prediction is not easy to perform. In contrast, the link-based recursive models are built based on the dynamic discrete choice framework RUST1987 and they are basically equivalent to discrete choice models based on the sets of all possible paths. Link-based recursive models have several advantages, making them attractive recently, for instance, the models can be consistently estimated and can be quickly and easily used for prediction. The first recursive model (i.e. the RL model) is proposed by FosgFrejKarl13. Several recursive models have been developed afterwards to deal with, for instance, correlations between path utilities MaiFosFre15,Mai_RNMEV,MaiBasFre15_DeC, dynamic networks de2020RL_dynamic, stochastic time-dependent networks mai2021RL_STD, or discounted behavior oyama2017discounted, with applications in, for instance, traffic management BailComi08,Melo12 or network pricing zimmermann2021strategic. We refer the reader to zimmermann2020tutorial for a comprehensive review.
With the development of technology and the ease of accessing the Internet, onboard devices, e.g., GPS, are more commonly used for several purposes, including rote navigation for taxis, buses, trucks, passenger cars, utility vehicles, mobile phones, etc. From operations and management perspectives, GPS data can be utilized to understand route choice behaviors. In fact, GPS has been used for collecting data for route choice modelling from a very early time wolf1999accuracy. Some other early studies used GPS data for modelling the route choice of vehicles gilani2005automatically, cyclists and pedestrians krizek2007detailed,dill2008understanding,menghini2010route,broach2012cyclists, transit oliveira2010improving, and commuting papinski2009exploring. While collecting data by GPS brings many advantages, it is frequent that GPS data is missing or incomplete, leading to difficulties for analysis wang2011challenges,shen2014review. Consequently, researchers have developed various approaches to overcome this issue. For example, smith2003exploring introduce a preliminary analysis of techniques for imputing missing data. The results revealed that the statistical technique outperformed the heuristics technique in terms of accuracy and performance. melville2004experiments used an ensemble method which is available under the DECORATE package and found DECORATE is more robust to missing data than the other two methods. Some researchers combined map-matching BierFrej08 or questionnaire surveys Ramos2015thesis which provide ground truth and improve the quality of the GPS data. We refer the reader to shen2014review for a review of GPS travel surveys, applications, and data processing. In the context of route choice modeling, to deal with the missing-data issue, to the best of our knowledge, existing studies only focus on data collection, i.e., improving GPS data and the map-matching, and we seem to be the first to take the estimation of route choice models with missing data into consideration.
Our work also closely relates to the EM framework dempster1977maximum,gentle1998algorithm,arcidiacono2011conditional, which has been popularly used for estimating statistical models with latent variables or unknown data observations. It should be noted that an EM algorithm does increase the likelihood value of the observed data after each iteration, but there is no guarantee that the sequence converges to a maximum likelihood estimator gentle1998algorithm. In addition, the method requires repeatedly solving MLE problems, which would be expensive, especially in the context of recursive route choice models since the estimation of these models involves solving several dynamic programs MaiFrejinger22.
Paper outline: This paper is organized as follows. Section (ref) presents the recursive models and the missing-data issue. Section (ref) presents our EM and Section (ref) presents the DC algorithms. Section (ref) reports numerical experiments, and finally, Section (ref) concludes the study and provides directions for future research.
Notations: Boldface characters represent matrices (or vectors), and $a_i$ denotes the $i$-th element of vector $\textbf{a}$, $A_{i,j}$ denotes the $(i,j)$-th element of matrix $\textbf{A}$, and $\textbf{A}_{:,j}$ denotes the $j$-th column of matrix $\textbf{A}$. We use $[m]$, for any $m\in \mathbb{N}_+$, to denote the set $\{1,\ldots,m\}$.
Consider a network ${\mathcal{G}} = ({\mathcal{A}},{\mathcal{V}})$, where ${\mathcal{A}}$ is the set of links and ${\mathcal{V}}$ is the set of nodes. Moreover, an absorbing state is associated with each destination by extending the network with dummy links $d$. The set of all links becomes $\widetilde{\mathcal{A}} = \mathcal{A} \cup \{d\}$. For each link $k\in {\mathcal{A}}$, let $A(k)$ be the set of outgoing links from $k$. For two links $k,a\in \widetilde{\mathcal{A}}$, the instantaneous utility $u(a|k) = v(a|k)+ \epsilon(a|k)$, where $v(a|k)$ is the deterministic part of the utility and $\epsilon(a)$ is a random term with zero mean. We assume that $v(a|k)$, for all $a\in A(k)$ to be negative for all links except the dummy link $d$.
For any link $k\in \mathcal{A}$, let $V(k)$ be the expected maximum utility (i.e., the value function) from the sink node of $k$ to the destination. In a recursive route choice model FosgFrejKarl13,MaiFosFre15, a choice of route is modelled as a sequence of link choices, in which at each link, the traveller selects an outgoing link by maximizing the sum of link utility and the expected maximum utility from that link to the destination. The selection of the next link can be written as in the following form $a = \text{argmax}_{a'\in A(k)} \left\{ v(a'|k)+ V(a') + \epsilon(a'|k)\right\}$. The value function $V(k)$ thus can be computed recursively as
Here, we note that the value function should depend on the destination $d$ but we omit the corresponding indices for notational simplicity. Different assumptions can be made for the random terms $\epsilon(a|k)$. If we assume that $\epsilon(a|k)$, $\forall a\in A(k)$, are i.i.d Extreme Value Type I with scale $\mu>0$ and independent of everything in the network, the recursive model becomes the recursive logit (RL) FosgFrejKarl13 in which the value function can be computed conveniently as
The above recursive formulation has an appealing property that it can be written as a system of linear equations. More precisely, if we define a vector $\textbf{Z}$ of size $|\widetilde{\mathcal{A}}|$ with elements $Z_k = \exp(V(k)/\mu)$, for all $k\in \widetilde{\mathcal{A}}$, a matrix ${\textbf{M}}$ of size $(|\widetilde{\mathcal{A}}|\times|\widetilde{\mathcal{A}}|)$ with elements $M_{ka} = \mathbb{I}(a\in A(k))\exp\left(v(a|k)/\mu\right)$, and $\textbf{b}$ is a vector of size $|\widetilde{\mathcal{A}}|$ with zero elements everywhere except $b_d=1$, then (ref) can be reformulated as a system of linear equations as $\textbf{Z} = {\textbf{M}}\textbf{Z} +\textbf{b}$. Such a system can be solved quickly using an iterative algorithm implemented in some existing linear solvers MaiFosFre15,MaiBasFre15_DeC.
It is well-known that the RL model retains the independence from irrelevance alternatives (IIA) property FosgFrejKarl13,MaiFosFre15, which would not hold in some real-life applications. To relax this property, MaiFosFre15 propose the NRL model based on the assumption that the random terms $\epsilon(a|k)$ have different scale parameters over links, i.e., $\epsilon(a|k) = \mu_k\epsilon'(a|k)$, $\forall k\in{\mathcal{A}}$, where $\epsilon'(a|k)$ are i.i.d extreme value type I for all $k,a\in\widetilde{\mathcal{A}}, a\in A(k)$. The value function in the NRL model can be defined recursively as
Then, by defining matrix ${\textbf{M}}$ of size $(|\widetilde{\mathcal{A}}|\times|\widetilde{\mathcal{A}}|)$ with entries $M_{ka} = \mathbb{I}(a\in A(k))\exp\left(v(a|k)/\mu_k\right)$, $\forall a,k\in\widetilde{\mathcal{A}}$, $\textbf{Z}$ of size $|\widetilde{\mathcal{A}}|$ with entries $Z_k = \exp(v(k)/\mu_k)$, the Bellman equation leads to the system of nonlinear equations $Z_k = \sum_{a\in A(k)} M_{ka}Z_a^{\phi_{ka}} +b_k$, $\forall k\in {\mathcal{A}}$. This system can be solved by a value iteration or least-squares approach MaiFrejinger22.
Recursive route choice models aforementioned can be estimated via MLE, which requires to compute log-likelihood values of some path observations. Given a complete path $\sigma = \{k_0,\ldots,k_{{\mathcal{T}}} = d\}$, the probability of choosing this path can be computed as \[ P(\sigma|{\pmb{\theta}}) = \prod_{t=0}^{{\mathcal{T}}-1} P(k_{t+1}|k_t;{\pmb{\theta}}), \] where ${\pmb{\theta}}$ is the model parameters to be estimated and $P(k_{t+1}|k_t;{\pmb{\theta}})$ is the probability of moving from link $k_t$ to $k_{t+1}$. The link choice probabilities $P(a|k)$, $\forall a\in A(k)$ can be further computed using the value function $V(k)$, $\forall k\in {\mathcal{A}}$. For instance, for the RL model, we can compute $P(a|k)$ as \[ P(a|k) = \exp\left(\frac{1}{\mu}\left(v(a|k) + V(a) - V(k)\right)\right), \] and for the NRL, the link choice probabilities can be computed as \[ P(a|k) = \exp\left(\frac{1}{\mu_k}\left(v(a|k) + V(a) - V(k)\right)\right). \] In the case that path observation $\sigma$ is not complete, i.e., there is at least one pair of successive link observations $(k_t,k_{t+1})$ such that $k_{t+t} \notin A(k_t)$, the computation of $P(\sigma)$ is much more challenging. To elaborate this point, let consider a link choice probability $P(a|k)$ where link $a$ is not connected to $k$, i.e., $a\notin A(k)$. Basically, $P(a|k)$ is referred to the probability of reaching link $a$ from $k$ conditional on the current parameters and the destination $d$. This probability can be written as
where $\Omega^{ka}$ is the set of all possible paths going from $k$ to $a$ and $P(\sigma)$ is the probability of (complete) path $\sigma$, which can be computed through link choice probabilities of connected links. Enumerating all paths from $k$ to $a$ is typically not feasible, in particular with real networks. So, a direct approach to deal with missing data is practical. In the next sections, we present our ways to address this issue, i.e., one is based on the EM scheme dempster1977maximum, and the other one is a direct approach to compute the log-likelihood function of incomplete observations.
We explore the EM method dempster1977maximum,arcidiacono2011conditional to handle the data missing issue. Such an EM algorithm will alternate an expectation (E) step, which creates a function for the expectation of the log-likelihood (LL) evaluated using the current estimate for the parameters ${\pmb{\theta}}$, and a maximization (M) step that maximizes the expectation function created by the E step. In the following, we show how this method can be applied to estimate recursive route choice models with incomplete observations.
We first denote by $\Phi = \{\sigma_1,\ldots,\sigma_N\}$ the set of $N$ observed paths. This set may contain some incomplete paths. We further denote by ${\pmb{\Gamma}} = \{\Gamma_1,\ldots,\Gamma_N\}$ the random variables representing the hidden parts of the observations. Each $\Gamma_n$ corresponds to $\sigma_n$ such that if $\gamma_n$ is a realization of $\Gamma_n$, then $(\sigma_n,\gamma_n)$ forms a complete path. Given a realization ${\pmb{\gamma}} = \{\gamma_1,\ldots,\gamma_N\}$ of ${\pmb{\Gamma}}$, we define $L({\pmb{\theta}}; \Phi,{\pmb{\gamma}})$ be the LL function of the complete observations $\Phi \cup {\pmb{\gamma}}$. To perform the two EM steps, we define $R({\pmb{\theta}}|\bar{{\pmb{\theta}}})$ as the expected value of the LL function of the observations w.r.t the current conditional distribution of ${\pmb{\Gamma}}$ given the observed set $\Phi$ and the current estimate $\bar{{\pmb{\theta}}}$ of the model parameters ${\pmb{\theta}}$. The function $R({\pmb{\theta}}|\bar{{\pmb{\theta}}})$ has the following form
Since $(\sigma_n,\gamma_n)$ forms a complete path for all $n\in [N]$, $P( \sigma_n,\gamma_n|{\pmb{\theta}})$ can be written as a product of choice probabilities of connected links. To compute the expectation in (ref), we consider $\sigma_n$ as a set of pairs of links and we separate each $\sigma_n$ into two sets of pairs of links $\sigma^1_n$ and $\sigma^2_n$ such that $\sigma^1_n$ contains all the connected pairs of links, i.e., $\sigma^1_n = \{(g,h)\in\sigma_n|\ h\in A(g)\} $ and $\sigma^2_n = \sigma_n\backslash \sigma^1_n$. This allows us to write (ref) as
where $\Gamma^{n,\varsigma}$ is a random variable representing a complete path connecting the pair of unconnected pair $\varsigma\in \sigma^2_n$, $ {\mathbb{E}}_{\Gamma^{n,\varsigma}|\varsigma,\bar{{\pmb{\theta}}}}\big[\ln P( \Gamma^{n,\varsigma}|{\pmb{\theta}})\big]$ is an expectation over $\Gamma^{n,\varsigma}$, and $P(\Gamma^{n,\varsigma}|{\pmb{\theta}})$ is the probability of complete path $\Gamma^{n,\varsigma}$. Since $\Gamma^{n,\varsigma}$ is complete, $P(\Gamma^{n,\varsigma}|{\pmb{\theta}})$ can be written as the product of choice probabilities of connected links.
An EM algorithm will perform alternatively E and M steps until converging to a fixed-point solution. At an E step and under a current estimate ${\pmb{\theta}}^t$, we need to create function $R({\pmb{\theta}}|{\pmb{\theta}}^t)$ based on $\Phi$ and the current parameters ${\pmb{\theta}}^t$. Even though there may be infinite numbers of complete paths $\Gamma$ connecting a pair of links $\varsigma \in \sigma^2_n$, we can approximate the expectation by sampling as
where $\gamma^{n,\varsigma}_1,\ldots,\gamma^{n,\varsigma}_S$ are $S$ realizations of $\Gamma^{n,\varsigma}$, $P(\gamma^{n,\varsigma}_s|\varsigma,{\pmb{\theta}}^t)$ is the probability of the path $\gamma^{n,\varsigma}_s$ given a pair of links $\varsigma$ and conditional on the current estimate ${\pmb{\theta}}^t$, and ${\mathcal{P}}^{n,\varsigma}({\pmb{\theta}}^t)=\sum_{s\in [S]} P(\gamma^{n,\varsigma}_s|\varsigma,{\pmb{\theta}}^t)$. Here, for ease of notation, we assume that the numbers of samples $S$ are the same over $\varsigma$. If it is not the case, we always can generate some dummy paths with zero probabilities. Typically, when the sample size $S$ is large enough, the left-hand side of (ref) will converge to the true $R({\pmb{\theta}}|{\pmb{\theta}}^t)$ with probability one due to the law of large number.
To sample from the conditional distribution of $\Gamma^{n,\varsigma}$ given a pair $\varsigma = (k,a)$, $a\notin A(k)$, and parameters ${\pmb{\theta}}^t$, we can generate paths that connect these two links according to the link probabilities $P(h|g;{\pmb{\theta}}^t)$, $\forall h\in A(g)$. The probability $P(\gamma^{n,\varsigma}|\;{\pmb{\theta}}^t)$ can be computed by multiplying the link choice probabilities $P(h|g;{\pmb{\theta}}^t)$ for every connected pair $(g,h)$ on $\gamma^{n,\varsigma}$. In summary, we describe the EM algorithm in Algorithm (ref) below.
The EM algorithm iterative performs the E and M steps and stops when it converges to a fixed point solution, i.e., ${\pmb{\theta}}^t$ is sufficiently closed to ${\pmb{\theta}}^{t+1}$. It generates a sequence $\{{\pmb{\theta}}^1,{\pmb{\theta}}^2,\ldots\}$ such that the marginal likelihood of the data $L(\Phi|{\pmb{\theta}}) = \int L({\pmb{\theta}};\Phi,\Gamma)d\Gamma$ is non-decreasing. In general, the algorithm does not guarantee that the sequence $\{{\pmb{\theta}}^t\}$ converges to a maximum likelihood estimator, but typically, the derivative of the likelihood function will be arbitrarily close to zero at the fixed point dempster1977maximum,arcidiacono2011conditional.
The M step of the EM algorithm involves solving a continuous maximization problem $\max_{{\pmb{\theta}}}\hat{R}({\pmb{\theta}}|{\pmb{\theta}}^t)$. The computation of $\hat{R}({\pmb{\theta}}|{\pmb{\theta}}^t)$ and its gradient w.r.t. ${\pmb{\theta}}$ can be done by computing the value function $V$ and its Jacobian. As mentioned earlier, these value functions can be obtained via solving systems of linear equations for the RL model, and via value iteration for the NRL model. In the context of RL, it can be shown that $\widehat{R}({\pmb{\theta}}|{\pmb{\theta}}^t)$ can be written as a sum of the logarithm of the choice probabilities of several complete paths. As a result, if the link utilities $v(a|k)$, $\forall a,k\in\widetilde{\mathcal{A}}$ are linear in ${\pmb{\theta}}$, then $\widehat{R}({\pmb{\theta}}|{\pmb{\theta}}^t)$ is concave in ${\pmb{\theta}}$, making the maximization problem in Algorithm (ref) tractable, noting that this linear-in-parameters setting is common in most of the existing route choice modeling studies Prato2012,FosgFrejKarl13,MaiFosFre15. We state this result in the following proposition.
\proof For notational simplicity, let us define $p^{n,\varsigma}_s = P(\gamma^{n,\varsigma}_s|\varsigma,{\pmb{\theta}}^t)/{\mathcal{P}}^{n,\varsigma}({\pmb{\theta}}^t)$, for all $n \in [N]$, $s\in [S]$, and $\varsigma\in \sigma^2_n$. For any $n,\varsigma$ we have $\sum_{s\in [S]} p^{n,\varsigma}_s = 1$. We now assume that each $\sigma^2_n$ contains $K_n$ disconnected pairs, denoted as $\varsigma^n_1,\ldots,\varsigma^n_{K_n}$ and let $${\mathcal{S}}^n = \underbrace{[S]\times \ldots\times [S]}_{K_n\text{ times}},$$ recalling that $[S] = \{1,2,\ldots,S\}$. We now write $\widehat{R}({\pmb{\theta}}|{\pmb{\theta}}^t)$ as
We now see that, for any $\varsigma^n_i\in \sigma^2_n$, $i\in [K_n]$, and $s\in [S]$, the term $\ln P( \gamma^{n,\varsigma^n_i}_s|{\pmb{\theta}}) p^{n,\varsigma^n_i}_s$ as
Moreover, we have
We now can combine (ref)-(ref) to have
where $P\left((\sigma_n,\gamma^{n,\varsigma^n_1}_{s_1},...,\gamma^{n,\varsigma^n_{K_n}}_{s_{K_n}} )|{\pmb{\theta}}\right)$ is the probability of the (complete) path $(\sigma_n,\gamma^{n,\varsigma^n_1}_{s_1},...,\gamma^{n,\varsigma^n_{K_n}}_{s_{K_n}} )$. FosgFrejKarl13 show that the choice probability of a complete path can be written as \[ P\left((\sigma_n,\gamma^{n,\varsigma^n_1}_{s_1},...,\gamma^{n,\varsigma^n_{K_n}}_{s_{K_n}} )|{\pmb{\theta}}\right) = \frac{\exp(v(\sigma_n,\gamma^{n,\varsigma^n_1}_{s_1},...,\gamma^{n,\varsigma^n_{K_n}}_{s_{K_n}} |{\pmb{\theta}}))}{\sum_{\tau \in \Omega^n} \exp(v(\tau|{\pmb{\theta}}))}, \] where $v(\tau|{\pmb{\theta}})$ is the deterministic utility of path $\tau$, i.e., $v(\tau) = \sum_{(k,a)\in\tau} v(a|k)$ and $\Omega^n$ is the set of all possible paths that share the same destination with observation $n${-th}. Under the assumption that $v(a|k)$ is linear in ${\pmb{\theta}}$, $v(\tau)$ is also linear in ${\pmb{\theta}}$. As a result, we can write
The first term of (ref) is linear in ${\pmb{\theta}}$ and the second term has a log-sum-exp convex form of a geometric program, thus is convex in ${\pmb{\theta}}$. So, $\ln P\left((\sigma_n,\gamma^{n,\varsigma^n_1}_{s_1},...,\gamma^{n,\varsigma^n_{K_n}}_{s_{K_n}} )|{\pmb{\theta}}\right)$ is concave, thus $\widehat{R}({\pmb{\theta}}|{\pmb{\theta}}^t)$ is concave, as desired. \endproof
If the choice model is the NRL, the LL and $\widehat{R}({\pmb{\theta}}|{\pmb{\theta}}^t)$ functions have highly nonlinear forms and are generally not concave. In this case, the M step of the EM algorithm may be difficult to handle globally. Moreover, it can be seen that the EM algorithm is generally expensive to perform, as it requires solving the maximization problem (ref) several times until getting a fixed point solution. Furthermore, each maximization problem (ref) requires sampling several paths between every unconnected pair of links to approximate the expectation ${R}({\pmb{\theta}}|{\pmb{\theta}}^t)$, which is also expensive.
In this section, we propose a practical way to exactly compute the likelihood of incomplete observations. We then present our decomposition-composition method to speed up the computation of the LL function as well as its gradients.
Given an unconnected pair $(k,a)$ such that the link observations between these two links are missing. To compute the LL of the path that contain pair $(k,a)$, we need to compute $P(a|k)$, i.e., the probability of reaching $a$ from $k$. As mentioned earlier, this can be done by enumerating all possible paths between these two links, which would be not possible with a real network. To do this in a tractable way, we define $\pi^a(s)$ as the probability of reaching link $a\in\widetilde{\mathcal{A}}$ from link $s\in{\mathcal{A}}$. The values of $\pi^a(s)$, $\forall s,a\in \widetilde{\mathcal{A}}$, satisfies the following recursive equations
So, if we define $\textbf{h}^a$ as a vector of size $|\widetilde{\mathcal{A}}|$ with all zero elements except $h^a_{a} = 1$ and a matrix $\textbf{Q}^a$ of size $(|\widetilde{\mathcal{A}}|\times|\widetilde{\mathcal{A}}|)$ with elements \[ Q^{a}_{s,s'} =
\] then we can write (ref) as
where $\pmb{\pi}^a$ is a vector of size $(|\widetilde{\mathcal{A}}|)$ with entries $\pi^a_s = \pi^a(s)$, $\forall s\in\widetilde{\mathcal{A}}$, and ${\textbf{I}}$ is the identity matrix of size $(|\widetilde{\mathcal{A}}|\times|\widetilde{\mathcal{A}}|)$. Proposition (ref) below shows that ${\textbf{I}}-\textbf{Q}^{a}$ is invertible, which guarantees the existence and uniqueness of the solutions to (ref).
So after computing $\pmb{\pi}^a$, we can obtain $P(a|k)$ as $P(a|k) = \pmb{\pi}^a_{k}$. We also need the gradients of $P(a|k)$ for the maximum likelihood estimation. Such gradient values can be obtained through taking the Jacobian of $\pmb{\pi}^a$ w.r.t. parameter $\theta_j$ as $$ \frac{\partial \pmb{\pi}^a}{\partial \theta_j} = -({\textbf{I}}-\textbf{Q}^{a})^{-1}\frac{\partial \textbf{Q}^a}{\partial \theta_j}({\textbf{I}}-\textbf{Q}^{a})^{-1}\textbf{h}^a = -({\textbf{I}}-\textbf{Q}^{a})^{-1}\frac{\partial \textbf{Q}^a}{\partial \theta_j} \pmb{\pi}^a, $$ where $\partial \textbf{Q}^a/\partial \theta_j$ is a matrix of size $|\widetilde{\mathcal{A}}|\times |\widetilde{\mathcal{A}}|$ with entries $\partial Q^a_{s,s'}/\partial \theta_j$. This matrix can be obtained by taking gradients of $P(k'|k)$ for all $k',k\in \widetilde{\mathcal{A}},k'\in A(k)$. These further can be computed by taking the gradients of the value functions, and the gradients of the value function $V(k)$ can be computed by solving some system of linear equations FosgFrejKarl13,MaiFosFre15. So, the Jacobian of $\pmb{\pi}^a$ can be obtained by solving a series of systems of linear equations.
In the above section, we show that the LL of incomplete observations can be computed directly by solving one system of linear equations for each unconnected pair $(k_i,k_{i+1})$. This would be time-consuming if a path observation contains a large number of such unconnected pairs. We show in the following that it is possible to obtain the probabilities of all the unconnected pairs on paths that share the same destination by solving only one system of linear equations, instead of solving one linear system per each incomplete pair. The idea is to decompose each individual system of linear equations into two parts among which one part is the same over unconnected pairs. This allows us to compose all the individual linear systems into only one linear system.
For a path observation, assume that we observe $J$ incomplete pairs $\{(u_1,v_1)$,..., $(u_J,v_J)\}$, $(u_j,v_j) \in \widetilde{\mathcal{A}}\times \widetilde{\mathcal{A}}$, $\forall j=1,\ldots,J$. We define a matrix $\textbf{Q}^0$ of size $((|\widetilde{\mathcal{A}}|+1)\times (|\widetilde{\mathcal{A}}|+1)$ with entries \[
\] That is, the last row and last column of $\textbf{Q}^0$ are all-zero vectors. {We also define a matrix $\textbf{H}$ of size $(|\widetilde{\mathcal{A}}|+1)\times K$ with entries} \[ H_{s,k} =
\] The following theorem show that one can obtain all the probabilities $P(v_j|u_j)$, $j=1,\ldots,J$ by solving only one system of linear equations. Let we first denote by $\widetilde{{\textbf{I}}}$ as an identity matrix of size $(|\widetilde{\mathcal{A}}|+1)\times (|\widetilde{\mathcal{A}}|+1)$.
\proof To prove the invertibility, Let $\widehat{\textbf{Q}} = ({\textbf{Q}}^0)^{\mbox{\tiny T}}$, we can see that $\widehat{Q}_{s,s'}\geq 0 $ for all $s,s'\in \widetilde{\mathcal{A}}$ and $ \sum_{s'\in A(s)} \widehat{Q}_{s,s'} = 1, \;\forall s\in \widetilde{\mathcal{A}} $ and $\sum_{s'\in A(s)} \widehat{Q}_{s,s'}=0 $ if $s = |\widetilde{\mathcal{A}}|+1$. Thus $\widehat{\textbf{Q}}$ is a sub-stochastic matrix. Moreover, $\widehat{\textbf{Q}}$ contains no recurrent class, thus $(\widetilde{{\textbf{I}}} - \widehat{\textbf{Q}})$ is invertible. Moreover $(\widetilde{{\textbf{I}}} - \widehat{\textbf{Q}})^{\mbox{\tiny T}} = (\widetilde{{\textbf{I}}} - {\textbf{Q}}^0)$, thus $\widetilde{{\textbf{I}}} - \textbf{Q}^0$ is invertible, as desired.
For the second claim of the theorem, let us consider a pair of missing data $(u_j,v_j)$. We create an artificial link $r$ such that it is impossible to reach $r$ from any other link in $\widetilde{\mathcal{A}}$ and from $r$ we can only reach $u_k$ with probability 1. Let $\pi^j(s)$, for all $s\in \mathcal{A}$ be the probability of reaching $s$ from $u_j$, i.e., $P(s|u_j)$. We see that $\pi^j(s)$, $s\in \mathcal{A}$, satisfy the following recursive equations.
Now, we define a matrix $\textbf{Q}^j$ of size $(|\widetilde{\mathcal{A}}|+1) \times (|\widetilde{\mathcal{A}}|+1)$ with elements $\textbf{Q}^j_{s,s'} = P(s|s')$, $\forall s,s'\in \widetilde{\mathcal{A}}\cup\{r\}, s\in A(s')$, and a vector $\textbf{q}$ of size $|\widetilde{\mathcal{A}}|+1$ with all zero elements except $q_r = 1$. We also number $r$ as $|\widetilde{\mathcal{A}}|+1$, so the last column and row of $\textbf{Q}^j$ and the last element of $\textbf{q}$ correspond to link $r$. (ref) can be written as
Note that the last row of $\textbf{Q}^j$ is an all-zero vector as $r$ is not reachable by any other links, and the last column of $\textbf{Q}^j$ has all zero entries except the element that corresponds to link $u_j$. We decompose (ref) as \[
\] where for (a) we separate $\textbf{Q}^j$ into two parts; the first part contains all the columns of $\textbf{Q}^j$ except the last column and the second part only contains the last column of $\textbf{Q}^j$. We now note that the last column of $\textbf{Q}^j$ have zero element except the one corresponding to $u_j$. Moreover $\pi^j_{|\widetilde{\mathcal{A}}|+1} = 1$, thus $(Q^j_{1,|\widetilde{\mathcal{A}}|+1},\ldots, Q^j_{|\widetilde{\mathcal{A}}|+1,|\widetilde{\mathcal{A}}|+1})^{\mbox{\tiny T}} \pi^j_{|\widetilde{\mathcal{A}}|+1}$ is a vector of zero elements except the one at $u_k$, which is equal to 1. Thus, \[
\pi^j_{|\widetilde{\mathcal{A}}|+1} + q \] is also the $j$-th column of matrix $\textbf{H}$ defined above. So, we have \[ \textbf{Q}^j\pmb{\pi}^j + \textbf{q} = \textbf{Q}^0\pmb{\pi}^j + \textbf{H}_{:,j}, \] where $\textbf{H}_{:,k}$ is the $j$-th column of $\textbf{H}$. So, $\pmb{\pi}^j$ is a solution to the following system \[ (\widetilde{{\textbf{I}}}-\textbf{Q}^0)\pmb{\pi}^j = \textbf{D}_{:,j}, \] where $\widetilde{{\textbf{I}}}$ is an identity matrix of size $(|\widetilde{\mathcal{A}}|+1)\times (|\widetilde{\mathcal{A}}|+1)$. This also means that if matrix $\pmb{\Pi}$ is a solution to the system of linear equations $(\widetilde{{\textbf{I}}}-\textbf{Q}^0)\pmb{\Pi} = \textbf{D}$ or $\pmb{\Pi} = (\widetilde{{\textbf{I}}}-\textbf{Q}^0)^{-1}\textbf{H}$, then \[ \pmb{\pi}^j = \pmb{\Pi}_{:,j}\text{ and } {\Pi}_{v_j,j} = \pmb{\pi}^j(v_j) = P(v_j|u_j), \] which is the desired result. \endproof
For the MLE, the gradients of $P(v_j|u_j)$ are required. The gradient of $P(v_j|u_j)$ w.r.t. a parameter $\theta_t$ can also be obtained by computing the Jacobian of $\pmb{\Pi}$, which can be obtained by solving the following linear system
Since $\widetilde{{\textbf{I}}}-\textbf{Q}^0$ is invertible, the linear system in (ref) always has a unique solution.
Algorithm (ref) below describes the basic steps of our Decomposition-Composition (DC) algorithm (the term “Decomposition-Composition” refers to the fact that we decompose matrix $\textbf{Q}^j$ into two parts where one of them is the same over unconnected pairs $(u_j,v_j)$, and then compose all the individual linear systems into only one system of linear equations).
The total number of linear systems to be solved in Algorithm (ref) to compute the probabilities of unconnected links is $(T+1)N^{\textsc{dest}}$, where $N^{\textsc{dest}}$ is the number of destinations and $T$ is the size of ${\pmb{\theta}}$. Here we note that if the model specification involves an origin-destination specific attribute, e.g., the link size attribute FosgFrejKarl13, then the total number of linear systems to be solved depends on the number of origin-destination pairs $N^{\textsc{od}}$ in the observations, i.e., $(T+1)N^{\textsc{od}}$. Clearly, the number of linear systems to be solved does not depend on the number of incomplete pairs of links in the observations. Moreover, Algorithm (ref) can be implemented in a parallel manner, which would help to speed up the computation. If the data is complete (no missing segment), then we just need to remove Step (iii) from the algorithm.
It should be noted that Step (iii) of Algorithm (ref) requires to solve some systems of linear equations that all involve the matrix ${\textbf{I}}-\textbf{Q}^0$. The $\textbf{L}\textbf{U}$ factorization is also a convenient approach to achieve good performance. Technically speaking, one can firstly decompose ${\textbf{I}}-\textbf{Q}^0$ into a lower triangular matrix $\textbf{L}$ and an upper triangular matrix $\textbf{U}$ and use these two matrices to solve the corresponding linear system.
To perform Step (iii) of Algorithm (ref), we need to compute $\partial P(a|k)/\partial \theta_i$ for any parameter $\theta_i$, for any $a\in A(k)$. $P(a|k)$ can be computed using the utility $v(a|k)$ and value function $V(a)$ and $V(k)$. For instance, for the RL model, $P_{ka} = \exp\left(\frac{1}{\mu}(v(a|k)+V(a)-V(k))\right)$. By taking the logarithm of $P_{ka}$ and taking the derivatives on the both sides, the derivative of $P(a|k)$ w.r.t a parameter $\theta_i$ can be easily obtained, for the RL model FosgFrejKarl13, as \[ \frac{\partial P_{ka}}{\partial \theta_i} = \frac{P_{ka}}{\mu}\left( \frac{\partial v(a|k)}{\partial \theta_i}+\frac{\partial V(a)}{\partial \theta_i} - \frac{\partial V(k)}{\partial \theta_i} \right), \] For the NRL model MaiFosFre15, the derivatives have more complicated formulas, as we assume that the scales $\mu_k$ are also functions of the model parameters ${\pmb{\theta}}$. \[ \frac{\partial P_{ka}}{\partial \theta_i} = P_{ka}\left(\frac{1}{\mu_k} \left(\frac{\partial v(a|k)}{\partial \theta_i}+\frac{\partial V(a)}{\partial \theta_i} - \frac{\partial V(k)}{\partial \theta_i} \right) - \frac{1}{\mu_k^2}\frac{\partial \mu_k}{\partial \theta_i} \left(v(a|k)+V(a)-V(k)\right)\right). \] Moreover, as shown in previous work MaiFosFre15,FosgFrejKarl13, the Jacobian of $V$ can be computed conveniently by solving a system of linear equations.
As discussed above, we have argued that the EM algorithm would be expensive as it requires several MLE problems until getting a fixed point solution. On the other hand, the DC algorithm requires to solve only one MLE problem, but at each iteration, we need to solve several systems of linear equations (or dynamic programs) to obtain the probabilities of unconnected pairs of links. So, it is to be expected that each iteration of the DC algorithm would be less expensive than an iteration of the EM, but more expensive than an iteration of the classical nested fixed-point algorithm used in previous work for the case of complete data MaiFrejinger22.
In this section, we present our numerical results with the RL FosgFrejKarl13 and the NRL models MaiFosFre15. We use a dataset collected in Borl\"ange, Sweden, which has been used in several prior route choice studies FosgFrejKarl13,MaiBasFre15_DeC,MaiFosFre15. The trip observations in the data are complete (there is no unconnected link observation) and to evaluate how the proposed algorithms handle the issue of data incompleteness, we will randomly remove links from the full observations.
We describe the dataset and our experimental settings in the following. The dataset is based on GPS observations of car trajectories collected in Borl\"ange, Sweden (a network composed of 3,077 nodes and 7,459 links). The sample consists of 1,832 trips corresponding to simple paths with a minimum of five links. Moreover, there are 466 destinations, 1,420 different origin-destination (OD) pairs and more than 37,000 link choices in this sample. We use the RL model specification of FosgFrejKarl13 and the NRL model specification of MaiFosFre15. The same four attributes as in FosgFrejKarl13,MaiFosFre15 are used for the instantaneous utilities, i.e., link travel time $TT(a)$ of link $a$, left turn dummy $LT(a|k)$ that equals one if the turn angle from $k$ to $a$ is larger than 40 and less than 177, U-turn dummy $UT(a|k)$ that equals one if the turn angle is larger than 177, and link constant $LC(a)$ that is equal to 1 for every link $a\in \mathcal{A}$. For the scales $\mu_k$ of the NRL model, we also use three attributes as in MaiFosFre15 for the scale specification $\mu_k$, which are travel time (TT), Link Size (LS) and the number of outgoing links (OL).
We generate datasets of missing information by removing links from the complete trip observations. To evaluate the performance of our proposed algorithms (i.e. the EM and DC algorithms), we will run them with the generated datasets and use the parameter estimates obtained to compute LL values based on the set of complete observations. We will compare the EM and DC algorithms with two baselines based on the standard estimation algorithm used in previous work, i.e. the nested-fixed point (NFXP) algorithm RUST1987. The first baseline is based on the complete dataset and the second one is only based on the connected links in the generated datasets. More precisely, for the second baseline, given an incomplete trip observation $\sigma = \{k_1,\ldots,k_L\}$, we compute the LL function based on $P(k_{i+1}|k_i)$ if $k_i$ and $k_{i+1}$ are connected, i.e., $k_{i+1} \in A(k_i)$, for $i\in [L-1]$. The difference between the two proposed algorithms (EM and DC) and the second baseline is that the EM and DC algorithms take into consideration the probabilities of unconnected links, whereas the baselines ignore them. As a result, it is expected that the two baselines would run faster than the EM and DC algorithms. We denote the first baseline by NFXP-C (standing for the NFXP algorithm with complete data) and the second baseline by NFXP-I (standing for the NFXP algorithm with incomplete data).
We generate datasets of missing information by going through all the trip observations in the dataset and randomly removing links (except the origins and destinations) with a given probability $p$. That is, given a removing probability $p \in [0,1]$, we remove each link observation with a probability $p$. We vary $p$ from 0.1 to 0.9, where small $p$ values imply that the missing information is low, and large $p$ values mean that the missing information is high. For each removing probability $p$, 10 seeds are created. Thus, there are 10 incomplete datasets generated for each removing probability $p$. In total there are 90 datasets generated.
We estimate the RL and NRL models using the EM, DC, NFXP-I algorithms with each generated incomplete dataset, and the NFXP-C with the complete dataset (the one before link removals). We then use the parameter estimates to compute the LL values of the complete dataset. Since the objective is to maximize the LL of the complete observations, the higher the LL values, the better the algorithm is. This will allow us to evaluate how the proposed algorithms (EM and DC) and the baseline NFXP-I recover the missing information, as compared to the NFXP-C which works with the complete dataset and always returns the highest LL values.
For the EM algorithm, we select the convergence threshold as $\xi = 10^{-4}$ and the number of path samples for each unconnected pair of links as $S=5$. Here we note that a higher number of samples can be chosen but since the EM algorithm is expensive to execute, especially when the number of path samples is large (as we will show later), we restrict ourselves to $S=5$.
The algorithms are coded in Python using Pytorch library paszke2019pytorch. The experiments are conducted on a PC with processor Intel(R) Xeon(R) Silver 4116 CPU @ 2.10GHz with 64Gb RAM, equipped with a GPU NVIDIA Quadro P1000. We use the limited memory BFGS (L-BFGS) method Nocedal2006 for the MLE problems (i.e., searching over the parameter space of ${\pmb{\theta}}$ for the DC algorithm and performing the M step of the EM algorithm). This is a Quasi-Newton method that computes the Jacobian matrix updates with gradient evaluations, but only saves some last updates to save memory.
In this section, we present numerical comparisons of the proposed algorithms (EM and DC) and the two baselines (NFXP-I and NFXP-C), based on the RL and NRL models. We first compare the LL values given by the four approaches, with a note that NFXP-C is based on the complete data and is considered as ground truth, thus always gives the highest LL values. The higher LL values imply that the corresponding algorithm performs better in recovering missing information from incomplete datasets.
Figure (ref) reports the means and standard errors of the LL values given by the four approaches, where the solid curves represent the mean values, and the shaded areas represent the standard errors. We can see that, for both RL and NRL models, the DC always gives higher LL values, as compared to those given by the EM and NFXP-I. The EM seems to be better than NFXP-I in some cases, and worse in some other cases. For the RL model, all the approaches (DC, EM and NFXP-I) seem to perform well as $p$ increases. The gaps between the DC, EM and NFXP-I are small for $p\leq 0.6$ and become more significant when $p\geq 0.7$. This demonstrates the efficiency of the DC algorithm in handling highly incomplete data. For the NRL model, we however see that the DC, EM and NFXP-I seem to perform worse, as compared to the case of the RL model; the LL values drop more quickly as $p$ increases. This is to be expected as we know that the NRL model has more complicated forms and has more information to infer (i.e., the parameters of the link utilities and the scales $\mu_k$, $k\in \widetilde{\mathcal{A}}$). We however still see that the DC significantly outperforms the EM and NFXP-I approaches for $p\leq 0.7$ and is slightly better than EM and NFXP-I for $p\geq 0.8$. Moreover, EM is slightly better than NFXP-I for $p\leq 0.5$ and slightly worse than NFXP-I for $p=0.9$.
We also report the details of the means and standard errors of the LL values in Table (ref), where the values before “$\pm$” are the means and the values after “$\pm$” are the standard errors. The LL values reported for $p=0.0$ are also the LL values given by the NFXP-C (i.e., given by the NFXP with the complete data). It can be seen that the LL values given by the NRL model are significantly higher than those given by the RL model. This is to be expected as previous studies already show that the NRL model always performs better than the RL in terms of in- and out-of-sample fits. Moreover, we see that the standards errors given by the RL are much smaller than those given by the NRL model, for all the three algorithms, which can be explained by the fact that the NRL model yields a more complicated choice probability formulation and is more difficult to estimate. Furthermore, looking at the LL values across the three approaches EM, DC, and NFXP-I, we also see that DC does not only return larger mean LL values, but also yields small LL standard errors. This would be due to the fact that the DC approach computes the probabilities of unconnected segments exactly while the EM algorithm approximates them by sampling, and the NFXP-I simply ignores the missing segments. This observation demonstrates the efficiency and stability of the DC approach, as compared to the other ones.
In summary, our experiments show that DC denominates the other approaches (i.e., EM and NFXP-I) in recovering missing information, in the sense that the LL values given by the DC algorithm are always higher (thus are closer to the ground-truth values), as compared to those given by EM and NFXP-I. Moreover, our approaches seem to provide better results for the RL model than for the NRL model.
We now turn our attention to the computational costs required to run our algorithms, which is also a focus of this work. In Figure (ref) we plot the means (solid curves) and standard errors (shaded areas) of the running times of the three algorithm DC, EM, and NFXP-I, for the RL and NRL models, noting that the running times of NFXP-I and NFXP-C are generally similar. In each figure, the left-hand side reports the average running times per iteration and the right-hand side reports the average total running times. It can be seen that, per iteration, DC and EM are much faster than EM, and as expected, NFXP-I is faster than DC. Importantly, the running times (per iteration) of EM grow quickly as $p$ increases while the running times (per iteration) of DC and NFXP-I are stable. This is because when the removing probability $p$ increases, the dataset would contain more unconnected segments and the EM algorithm would need more samples to approximate the expected LL function, thus becoming more expensive. In contrast, the DC algorithm is developed in such a way that the number of systems of equations to be solved does not depend on the number of unconnected segments in the path observations. The number of linear systems to be solved also does not change as $p$ increases for the NFXP-I algorithm. Thus, an increase in $p$ would not make a significant impact on the running times of DC and NFXP-I.
In terms of total estimation time, it is interesting to see that the EM algorithm is faster than the DC for some small $p$ values (e.g., $p\leq 0.2$ for the RL and $p\leq 0.3$ for the NRL models), even though the running times per iteration of the EM are higher. This would be explained as follows. When $p$ is small, the number of unconnected segments in the observations is small, thus the EM algorithm only requires a small number of path samples and $\widehat{R}({\pmb{\theta}}|{\pmb{\theta}}^t)$ would give a good approximation to the expected LL function. Moreover, the form of the approximate expected LL function in (ref), even though looks complex, it shares a similar structure with the standard LL functions studied in previous works FosgFrejKarl13,MaiFosFre15. On the other hand, the LL function in the DC approach involves some exact probabilities of unconnected segments, thus would be more complicated. All these would be the reason for the fact that when $p$ is small, EM requires fewer iterations to converge and is faster, as compared to the DC algorithm.
When $p\geq 0.2$ for the RL model and $p\geq 0.3$ for the NRL model, the running times of EM start becoming higher than those of DC and increase fast. This is to be expected as the running times per iteration of EM also increases fast as $p$ increase. The total estimation times of the DC algorithm, as expected, are stable and seem not to be affected much by an increase in $p$. Moreover, looking across the two route choice models, we see that the total running times of DC are higher for the NRL model, as compared to those given by the RL model. This is to be expected as the NRL model is known to be more expensive to estimate than the RL model. However, it seems not to be the case for the EM approach. The standard errors of the running times seem to be higher for the NRL model, which is consistent with our above observations for the LL comparisons. We refer the reader to the appendix for the details of the running times.
In summary, in terms of running time, the DC and EM algorithms are more expensive to perform than the naive approach NFXP-I. The DC algorithm is generally much faster than the EM, especially when $p$ is large. More importantly, the running times of EM grow fast as the data-missing level $p$ increases, while an increase in $p$ does not make a significant impact on the running times of DC. All these indicate the practical tractability of the DC approach for handling highly missing data.
In this paper, we have studied the issue of incomplete observations in route choice modeling. we have developed two solution approaches, namely, the EM and DC algorithms, aiming to efficiently recover missing information from data. While the EM algorithm requires to solve several MLE problems and would be expensive to perform, the DC algorithm is based on the idea that we can exactly compute the probabilities of unconnected segments by solving systems of linear equations. The main advantage of our DC algorithm is that the number of system of linear equations to be solved to compute the probabilities of unconnected segments does not depend on the number of unconnected segments, thus this approach scales well when the number of unconnected segments increases. We provide numerical experiments based on a dataset collected from a real network, which showed that the DC algorithm outperforms the EM and other baselines in recovering missing information. Our numerical results also showed that the DC algorithm is generally faster than the EM approach, and the running times of the DC seem not to be affected by the amount of missing information in the data observations.
Our methods can be used with other recursive models developed in the literature, e.g., the dynamic RL de2020RL_dynamic or stochastic time-dependent RL mai2021RL_STD models, and may have applications beyond route choice modeling, for example, the estimation of activity-based models or general structural dynamic discrete choice models with missing data. These would shape some interesting directions for future work.
{This research is supported by Singapore Ministry of Education (MOE) Academic Research Fund (AcRF) Tier 1 grant (Grant No: 20-C220-SMU-010) to the first and second authors}