203: Exampville Destination Choice

Welcome to Exampville, the best simulated town in this here part of the internet!

Exampville is a demonstration provided with Larch that walks through some of the data and tools that a transportation planner might use when building a travel model.

import larch, numpy, pandas, os
from larch import P, X

Larch 5.2.11 /Users/jpn/Git/Larch/larch

In this example notebook, we will walk through the estimation of a tour destination choice model. First, let’s load the data files from our example.

hh, pp, tour, skims, emp = larch.example(200, ['hh', 'pp', 'tour', 'skims', 'emp'])
For this destination choice model, we’ll want to use the mode choice logsums we calculated previously from the mode choice estimation, but we’ll use these values as fixed input data instead of a modeled value.
We can load these logsums from the file in which they were saved. For this example, we can indentify that file using the larch.example function, which will automatically rebuild the file if it doesn’t exists. In typical applications, a user would generally just give the filename as a string and ensure manually that the file exists before loading it.
logsums_file = larch.example(202, output_file='/tmp/logsums.pkl.gz')
logsums = pandas.read_pickle(logsums_file)


We’ll replicate the pre-processing used in the mode choice estimation, to merge the household and person characteristics into the tours data, add the index values for the home TAZ’s, filter to include only work tours, and merge with the level of service skims. (If this pre-processing was computationally expensive, it would probably have been better to save the results to disk and reload them as needed, but for this model these commands will run almost instantaneously.)

raw = tour.merge(hh, on='HHID').merge(pp, on=('HHID', 'PERSONID'))
raw["HOMETAZi"] = raw["HOMETAZ"] - 1
raw["DTAZi"] = raw["DTAZ"] - 1
raw = raw[raw.TOURPURP == 1]
raw.index.name = 'CASE_ID'

The alternatives in the destinations model are much more regular than in the mode choice model, as every observation will have a similar set of alternatives and the utility function for each of those alternatives will share a common functional form. We’ll leverage this by using idca format to make data management simpler.

First, we’ll assemble some individual variables that we’ll want to use. We can build an array to represent the distance to each destination based on the "AUTO_DIST" matrix in the skims OMX file.

distance = pandas.DataFrame(
    data=skims.AUTO_DIST[:][raw["HOMETAZi"], :],

This command pulls the relevant row, identified by the "HOMETAZi" column in the raw data, into each row of a new DataFrame, which has a row for each case and a column for each alterative.

Note that the [:] inserted into the data argument is used to instruct the pytables module to load the entire matrix into memory, and then numpy indexing is used to actually select out the rows needed. This is a technical limitation of the pytables module and could theoretically be a very computationally expensive step if the skims matrix is huge relative to the number of rows in the raw DataFrame. However, in practice a single matrix from the skims file is generally not that large compared to the number of observations, and this step can be processed quite efficiently.

The logsums we previously loaded is in the same format as the distance, with a row for each case and a column for each alterative. To use the idca format, we’ll reshape these data, so each is a single column (i.e., a pandas.Series), with a two-level MultiIndex giving case and alternative respectively, and then assemble these columns into a single DataFrame. We can do the reshaping using the stack method, and we will make sure the resulting Series has an appropriate name using rename, before we combine them together using pandas.concat:

ca = pandas.concat([
], axis=1)
<class 'pandas.core.frame.DataFrame'>
MultiIndex: 198080 entries, (0, 1) to (15931, 40)
Data columns (total 2 columns):
distance    198080 non-null float64
logsum      198080 non-null float64
dtypes: float64(2)
memory usage: 3.6 MB

For our destination choice model, we’ll also want to use employment data. This data, as included in our example, has unique values only by alternative and not by caseid, so there are only 40 unique rows. (This kind of structure is common for destination choice models.)

<class 'pandas.core.frame.DataFrame'>
Int64Index: 40 entries, 1 to 40
Data columns (total 3 columns):
NONRETAIL_EMP    40 non-null int64
RETAIL_EMP       40 non-null int64
TOTAL_EMP        40 non-null int64
dtypes: int64(3)
memory usage: 1.2 KB

To make this work with the computational arrays required for Larch, we’ll need to join this to the other idca data. Doing so is simple, because the index of the emp DataFrame is the same as the alternative id level of the ca MultiIndex. You can see the names of the levels on the MultiIndex like this:

FrozenList(['CASE_ID', 'TAZ_ID'])

Knowing the name on the alternatives portion of the idca data lets us join the employment data like this:

ca = ca.join(emp, on='TAZ_ID')

Then we bundle the raw data along with this newly organized idca data, into the larch.DataFrames structure, which is used for estimation. This structure also identifies a vector of the alterative codes and optionally, names and choice identifiers. This structure can be attached to a model as its dataservice.

dfs = larch.DataFrames(
    alt_names=['TAZ{i}' for i in skims.TAZ_ID],
larch.DataFrames:  (not computation-ready)
  n_cases: 4952
  n_alts: 40
    - distance
    - logsum
    - TOURID
    - HHID
    - DTAZ
    - X
    - Y
    - INCOME
    - geometry
    - HHSIZE
    - HHIDX
    - AGE
    - WORKS
    - N_TOURS
    - HOMETAZi
    - DTAZi
  data_av: <populated>
  data_ch: DTAZ

Model Definition

m = larch.Model(dataservice=dfs)
m.title = "Exampville Work Tour Destination Choice v1"
m.quantity_ca = (
        + P.EmpRetail_HighInc * X('RETAIL_EMP * (INCOME>50000)')
        + P.EmpNonRetail_HighInc * X('NONRETAIL_EMP') * X("INCOME>50000")
        + P.EmpRetail_LowInc * X('RETAIL_EMP') * X("INCOME<=50000")
        + P.EmpNonRetail_LowInc * X('NONRETAIL_EMP') * X("INCOME<=50000")

m.quantity_scale = P.Theta

m.utility_ca = (
    + P.logsum * X.logsum
    + P.distance * X.distance

Model Estimation

req_data does not request {choice_ca,choice_co,choice_co_code} but choice is set and being provided
req_data does not request avail_ca or avail_co but it is set and being provided

Iteration 005 [Converged]

LL = -16366.83084793109

value initvalue nullvalue minimum maximum holdfast note best
EmpNonRetail_HighInc 0.903947 0.0 0.0 -inf inf 0 0.903947
EmpNonRetail_LowInc -0.997606 0.0 0.0 -inf inf 0 -0.997606
EmpRetail_HighInc 0.000000 0.0 0.0 0.000 0.0 1 0.000000
EmpRetail_LowInc 0.000000 0.0 0.0 0.000 0.0 1 0.000000
Theta 0.728573 1.0 1.0 0.001 1.0 0 0.728573
distance 0.006269 0.0 0.0 -inf inf 0 0.006269
logsum 1.200129 0.0 0.0 -inf inf 0 1.200129
EmpNonRetail_HighInc 0.903947
EmpNonRetail_LowInc -0.997606
EmpRetail_HighInc 0.000000
EmpRetail_LowInc 0.000000
Theta 0.728573
distance 0.006269
logsum 1.200129
stepsarray([1., 1., 1., 1., 1.])
message'Optimization terminated successfully.'

Model Visualization

For destination choice and similar type models, it might be beneficial to review the observed and modeled choices, and the relative distribution of these choices across different factors. For example, we would probably want to see the distribution of travel distance. The Model object includes a built-in method to create this kind of visualization.


The distribution_on_continuous_idca_variable has a variety of options, for example to control the number and range of the histogram bins:

m.distribution_on_continuous_idca_variable('distance', bins=40, range=(0,10))

Alternatively, the histogram style can be swapped out for a smoothed kernel density function:


Subsets of the observations can be pulled out, to observe the distribution conditional on other idco factors, like income.

    xlabel="Distance (miles)",
    header='Destination Distance, Very Low Income (<$10k) Households',

Destination Distance, Very Low Income (<$10k) Households

Save and Report Model

report = larch.Reporter(title=m.title)
report << '# Parameter Summary' << m.parameter_summary()

Parameter Summary

ParameterValueStd Errt StatNull Value
EmpRetail_HighInc0fixed value
EmpRetail_LowInc0fixed value
report << "# Estimation Statistics" << m.estimation_statistics()

Estimation Statistics

StatisticAggregatePer Case
Number of Cases4952
Log Likelihood at Convergence-16366.83-3.31
Log Likelihood at Null Parameters-18362.98-3.71
Rho Squared w.r.t. Null Parameters0.109
report << "# Utility Functions" << m.utility_functions()

Utility Functions

X.logsumThis is Data

X.distanceThis is Data

* log(
   + exp(
P.EmpRetail_HighIncexp(0) = 0
) *
X('RETAIL_EMP * (INCOME>50000)')This is Data

   + exp(
P.EmpNonRetail_HighIncexp(0.904) = 0.9039
) *
X('NONRETAIL_EMP*(INCOME>50000)')This is Data

   + exp(
P.EmpRetail_LowIncexp(0) = 0
) *
X('RETAIL_EMP*(INCOME<=50000)')This is Data

   + exp(
P.EmpNonRetail_LowIncexp(-0.998) = -0.9976
) *
X('NONRETAIL_EMP*(INCOME<=50000)')This is Data


The figures shown above can also be inserted directly into reports.

figure = m.distribution_on_continuous_idca_variable(
    xlabel="Distance (miles)",
    header='Destination Distance',
report << "# Visualization"
report << figure


Destination Distance