# Copula-based synthetic population generation

Pascal Jutras-Dubé<EMAIL_ADDRESS>Mohammad B. Al-Khasawneh
<EMAIL_ADDRESS>Zhichao Yang<EMAIL_ADDRESS>Javier Bas<EMAIL_ADDRESS>Fabian Bastin<EMAIL_ADDRESS>Cinzia Cirillo<EMAIL_ADDRESS>Department of Computer Sience and Operations Research, Université de Montréal,
CP 6128, Succ. Centre-Ville, Montreal, QC, H3C 3J7, Canada  CIRRELT,
Université de Montréal, CP 6128, Succ. Centre-Ville, Montreal, QC, H3C 3J7,
Canada  Department of Civil and Environmental Engineering, University of
Maryland, College Park, MD, 20742, USA  Department of Quantitative Economics,
Universidad Autónoma de Madrid, Madrid, 28049, Spain

###### Abstract

Population synthesis consists of generating synthetic but realistic
representations of a target population of micro-agents for the purpose of
behavioral modeling and simulation. We introduce a new framework based on
copulas to generate synthetic data for a target population of which only the
empirical marginal distributions are known by using a sample from another
population sharing similar marginal dependencies.  This makes it possible to
include a spatial component in the generation of population synthesis and to
combine various sources of information to obtain more realistic population
generators. Specifically, we normalize the data and treat them as realizations
of a given copula, and train a generative model on the normalized data before
injecting the information on the marginals. We compare the copulas framework
to IPF and to modern probabilistic approaches such as Bayesian networks,
variational auto-encoders and generative adversarial networks. We also
illustrate on American Community Survey data that the method proposed  allows
to study the structure of the data at different geographical levels in a way
that is robust to the peculiarities of the marginal distributions.

###### keywords:

Population synthesis , Generative modeling , Model transferability ,
Statistical learning algorithms , Copulas , Small areas.

##  1 Introduction

Population synthesis refers to models that aim at constructing artificial
datasets whose characteristics mimic that of a population of interest.
Population synthesizers can produce a set of agents with detailed
sociodemographic and socioeconomic information at both the individual and
household levels, maintaining the structural coherence of the population they
replicate. Therefore, these methodologies have been of great use wherever an
agent-based model was to be implemented, such as transportation models based
on micro-simulation [2, 49, 42, 28],  for which analyzing the spatial
implications of a given policy is capital.

In order to properly reflect the population properties it aims to  mimic, a
synthetic population must share the same joint distribution of the variables
that it encapsulates. In the context of microsimulation frameworks for travel
behavior, this property must hold for the geographic area of interest. Various
techniques have been proposed to achieve this goal. However, with a few
notable exceptions, as  in Barthelemy and Toint [6], most of them rely on a
sample of the target population (census data, travel survey), often costly to
obtain. As a result, the sample size can be limited, especially when working
at small geographical level ( i.e. census tract level). It is therefore
desirable to include a spatial component in the population synthesis
generation and to combine various information sources to obtain more realistic
population generators. However, such aspects of the problem have received very
limited attention in the transportation literature. Thus, we here introduce a
framework that can generate synthetic data for a target population of which
only the marginal or empirical marginal distributions are known at different
geographical levels ( i.e. state, county, Public Use Micro Areas [PUMA],
census tract) by using a sample from another population sharing similar
marginal dependencies. To achieve this objective, we combine the theory of
copulas with traditional machine learning (ML) generative modeling approaches
to separate the learning of the dependencies structure from the marginal
distributions, allowing its use with various populations that differ on the
marginals only. We compare several methods (Bayesian networks, variational
autoencoders, generative adversarial networks) alone or in combination with
copulas, and illustrate on American Community Survey (ACS) data that the
method proposed makes it possible to study the structure of the data at
different geographical levels in a way that is robust to the peculiarities of
the marginal distributions.

The rest of the paper is organized as follows. In Section 2, we review the
literature on population synthesis, with a focus on the different existing
methodologies recently proposed. Section 3 describes the methods applied in
this work, while the specifics of the experiments carried out, as well as the
results obtained, are detailed in Section 4. Finally, conclusions and future
research avenues are presented in Section 5.

##  2 Literature Review

The ultimate goal of a population synthesizer is to use in different processes
(modeling, optimization, simulation, etc.) the new information that it
generates from an available disaggregated dataset insufficient in size. There
are several notable examples in transportation where this exercise is
performed. Eluru et al. [21] develop a microsimulator called CEMUS that feeds
aggregated socioeconomic information into a population generator to produce
household and individual disaggregated synthetic datasets to model individual-
level activity-travel patterns. Bradley et al. [15] present a forecasting
model implemented by the Sacramento Area Council of Governments that relies
upon the simulation of residents’ full-day activity and travel schedule. Auld
et al. [3] use a dynamic simulator to feed the Transportation Analysis and
Simulation System (TRANSIMS) [57] in order to generate route choices. Ziemke
et al. [65] use a synthetic population to compute accessibility measures in
Nelson Mandela Bay in South Africa by using census and travel survey data.

A desirable property of any population synthesizer is to preserve the general
characteristics of the population, which leads to different technical
approaches. The most widely used, until recently, has probably been the
Iterative Proportional Fitting (IPF), initially proposed by Deming and Stephan
[19], and popularized in transportation by Duguay et al. [20]. IPF selects
households from the source sample trying to match some given marginals totals,
requiring a fitting stage and an allocation stage. In the fitting step, a
contingency table is computed from the seed table (the source sample) and the
marginals totals. In the allocation phase, households are randomly selected
from the seed table to match the frequency given in the contingency table.
There is an abundant literature on empirical applications of IPF —reviewed,
for instance, by Müller and Axhausen [42], as well as on its technical aspects
[1, 53, 4, 62, 28]. The method nevertheless presents several important flaws.
In this regard, Pritchard and Miller [49] address computational memory
restrictions, while Guo and Bhat [28] explore how to avoid sampling zero
issues. Ye et al. [62] consider simultaneously fitting different types of
agents, proposing a heuristic approach called Iterative Proportional Updating
(IPU) to overcome the disadvantages of the standard IPF. However, they fail to
accommodate the new synthetic information at multiple geographical resolutions
simultaneously, leading to a loss of representativeness. Konduri et al. [38]
extend their efforts, proposing an enhanced IPU algorithm that accounts for
constraints at different levels of spatial resolution when generating a
synthetic population.

However, as Farooq et al. [22] point out, fitting a contingency table to the
available data may entail errors if the information is not complete or has
been manipulated. In fact, these potential errors cannot be contrasted due to
absence of a complete real data set. Thus, these authors propose a Markov
Chain Monte Carlo (MCMC) simulation-based approach that uses partial views of
the joint distribution of the real population obtained from the census to draw
from it high-dimensional synthetic populations that would otherwise have been
impossible to produce using IPF.  MCMC technique was further explored by
Casati et al. [16] and Saadi et al. [51, 52].

Another category of methods for synthesizing populations, the Combinatorial
Optimization Methods, proceed differently. They divide an area into mutually
exclusive subareas where a distribution of the attributes of interest is
available. Then a sample taken over the whole population is fitted to the
given set of marginals for each subarea. Huang and Williamson [30] offer a
comparison of these two approaches in the context of the creation of small-
area microdata. However, the assumptions about the data that these two
families of methods require cannot always be met. This led Barthelemy and
Toint [6] to develop a new type of generator, applied to the case of Belgium.
Their technique works in a three-step process. First, a pool of individuals
pertaining to a certain household is generated; then, the household joint
distribution is estimated and stored in a contingency table; finally, the
synthetic households are constructed by randomly drawing individuals from the
pool of individuals, preserving the distribution computed in the second step.
They compare the results of this method with those of the extended IPF
described in Guo and Bhat [28], observing that the new generator performs
better. Another effort to offer an alternative approach is  due to Sun et al.
[59], who propose a hierarchical mixture model to generate representative
household structures in population synthesis. Their framework comprises a
probabilistic tensor factorization, a multilevel latent class model, and a
rejection sampling process. They test their procedure on the Household
Interview Travel Survey data of Singapore, being able to generalize the
associations among attributes as well as to reproduce structural relationships
among household members.

Yet, it is possible to observe a recent increasing trend in the substitution
of iterative fitting algorithms for the generation of synthetic populations by
alternative approaches. One of these new perspectives is the application of
Copula Generative Models. Copulas are a mathematical probability tool that
allows the modeling of random variables that have an intrinsic relationship to
each other. The dependence structure can be decoupled from the original
available information so  that any set of new synthetic data can be created
having that same structure. We refer to Genest and Favre [24] for a gentle
introduction, and to the books by Joe [33, 34] and Nelsen [43] for readers
seeking concrete mathematical proofs, theorems and derivations related to
copulas. The selection of a specific copula can remain challenging [36],
especially in the case of discrete copula.  Avramidis et al. [5] study the
NORTA correlation-matching problem for modeling dependence of discrete
multivariate distributions, while Nikoloulopoulos [44] reviews copula-based
approaches with discrete multivariate response data. Copulas have been
successfully applied in various domains, including finance and actuarial
sciences [17], transportation [7, 48, 50, 35, 11], water resources management
[10], and simulation modeling of arrival rates in call centers [46, 32].

Bayesian network (BN) is another alternative approach that has recently gained
momentum as a tool to perform population synthesis. This method encodes the
dependency relationships among predictors using a graphical model in which
nodes represent variables, links  capture their conditional dependencies, and
probability distributions are assigned to each node, conditional on its
parents. As expressed in the seminal work  by Sun and Erath [58], if the
conditional structure of the data generating process is known, inference can
be based on the underlying probabilities and the population can then be
generated accordingly by sampling from the joint distribution. Since this
information is typically not available, the authors propose that the graph
structure of the data may be learned through a scoring approach. More recent
examples of the use of BN for population synthesis are given by Zhang et al.
[64], who make use of traditional survey data and digital records of
networking and human behavior to generate connected synthetic populations
using BN; Hörl and Balac [31], from their part, generate a synthetic travel
demand based on open data.

On the other hand, other approaches to the generation of synthetic populations
fall more along the lines of ML. The most popular of these methodologies are
probably generative models [8], given their success in generating artificial
images. A generative model captures the dependency structure of the variables
of a dataset and trains a generator to reproduce samples preserving their
joint distribution. Among generative models, two families stand out:
variational autoencoders (VAE) [37] and generative adversarial networks (GAN)
[25]. VAE are unsupervised generative models composed of an encoder and a
decoder that are trained to minimize the reconstruction error between the
encoded/decoded data and the initial data. When the model is trained, it
learns the joint distribution of the data.  Borysov and Rich [12] use a
conditional VAE to study travel preference dynamics using a synthetic pseudo-
panel approach, while Boquet et al. [9] propose a VAE model to generate
traffic data. In particular, we use in this paper  a tabular VAE (TVAE) [60],
an adaptation of regular VAE to tabular data that applies the same
preprocessing but modifies the loss function.  On the other hand, GANs were
first introduced by Goodfellow et al. [25] as a framework in which two neural
networks, called generator (which produces the target output) and
discriminator (which distinguish true data from the output of the generator),
compete with one another to deceive and not be deceived. This process leads to
a model that captures the distribution and dependency of the features of a
dataset, and ultimately creates new synthetic information. Despite  their
great success in the computer vision and natural language processing fields,
there are only a few attempts of using GAN s in transportation. Yazdizadeh et
al. [61] use  a GAN to infer travel mode from GPS trajectory data and derive
trip information from travel survey. The work of Günthermann et al. [27]
delves along the same lines. On the contrary, Yin et al. [63] focus on the use
of deep generative models to generate synthetic travelers’ mobility patterns
that replicate the statistical properties of a sample of actual travelers.
Given the structure of our data,  in this paper, we use a particular GAN , as
in the case of VAE;, the so-called  conditional tabular GAN (CTGAN) [60]. It
was designed to adapt  a GAN to tabular data with a mix of potentially highly
imbalanced discrete columns and multimodal non-Gaussian-like continuous
columns.

##  3 Methodology

Consider a sample of size $N$ from a random vector $X$, from which we want to
generate a synthetic population of any arbitrary size. As  highlighted in the
literature review, many techniques can be used to generate a synthetic
population , yet we will  only select, in addition to IPF, the most promising
approaches.

We first briefly introduce Bayesian networks, conditional tabular generative
adversarial network (CTGAN), and tabular variational autoencoders (TVAE) [60],
that will serve as our benchmark methods. We then discuss copulas and their
link to normalization procedures in ML, and conclude with their use  in the
generation of synthetic populations.

###  3.1 Bayesian networks

A Bayesian network (BN) is a probabilistic model that represents  a set of
random variables $X_{i}$, $i=1,\ldots,d$ , and their conditional dependencies
in the form of a directed acyclic graph $G$. Each node corresponds to one
random variable, and the arcs express the conditional dependencies between
them. If no path exists between two nodes, the corresponding random variables
are independent. The joint probability distribution of $X$ can be decomposed
as the product of the marginal distributions of $X_{i}$, conditioned on its
parents $\Pi_{i}$, $i=1,\ldots,d$, using the chain rule

$P(X\mid G,\Theta)=\prod_{i=1}^{d}P(X_{i}\mid\Pi_{i},\theta_{i}),$

where $\theta_{i}$ are the parameters of the distribution of $X_{i}$ and
$\Theta=\left\\{\theta_{1},\ldots,\theta_{d}\right\\}$.

Both the graph structure $G$ and the parameters $\Theta$ are learned. This
process is often referred to as  structural learning and has two stages: model
selection and model optimization. In the selection stage, a score function is
used to quantify how well a hypothetical graph structure fits the data. The
minimum description length (MDL) is a popular scoring function consisting of
two components that estimate the structural complexity and the likelihood of
the data given the model [39]. For discrete BN, the variables conditioned on
their parents are assumed to follow multinomial distributions and their
parameters are estimated by maximum likelihood estimation. The goal of the
model optimization stage is to identify the hypothetical structure with the
highest score. However, the evaluation of all possible graph structures has
unrealistic time complexity [18]. Instead, we can apply a heuristic search
that greedily chooses a topological ordering of the variables, and optimally
identifies the best parents for each variable given this ordering [29], as
illustrated in Figure 1, depicting a Bayesian network built on our dataset.

Figure 1: Example of Bayesian Network found by the greedy search algorithm.

Once the structural learning is completed, a synthetic population can be drawn
from the the factorized joint probability distribution defined by the Bayesian
network. More precisely, we generate values via rejection sampling,  and refer
the reader to Sun and Erath [58] for more details on the application of BN in
population synthesis.

###  3.2 Conditional tabular generative advers arial network

Generative adversarial networks (GANs) [25] are a class of unsupervised ML
algorithms in which a generative model is opposed to a discriminative model
that learns to determine whether a sample comes from the model distribution or
from the data distribution. The competition causes both parties to improve
until the generative model can generate false data that is indistinguishable
from the real data.

The generator $G_{\theta}$ is a neural network with parameters $\theta$ that
approximates the distribution of data $p_{X}$ by generating synthetic samples
$G_{\theta}(Z)$ from a prior random distribution $Z\sim p_{Z}$. The
discriminator $D_{\phi}$ is a neural network with parameters $\phi$ whose
scalar output $D_{\phi}(x)$ represents the probability that an example $x$
comes from the real data rather than from the distribution of the generator.
$D_{\phi}$ is trained to maximize the probability of assigning the correct
class whether an example was sampled from the training set or from
$G_{\theta}$. $G_{\theta}$ is simultaneously trained to minimize the
probability of $D_{\phi}$ assigning the correct class to its output

$\min_{\theta}\max_{\phi}{\mathbb{E}}_{X\sim p_{X}}[\log
D_{\phi}(X)]+{\mathbb{E}}_{Z\sim p_{Z}}[\log(1-D_{\phi}(G_{\theta}(Z)))].$

GANs have achieved impressive performance in generating high-quality synthetic
images but were not explicitly designed to manipulate tabular data. This led
to the development of conditional tabular GANs (CTGAN) [60] which are capable
of handling tabular data with a mix ture of potentially highly imbalanced
discrete columns and multimodal  non-Normal continuous columns.

A CTGAN uses mode-specific normalization to deal with columns with complicated
distributions, which is often the case for continuous features in tabular
data. For a continuous column, mode-specific normalization first estimates the
number of modes with a variational Gaussian mixture model (VGM). Then, for
each value in the column, it calculates the probability that it comes from
each mode with the distributions given by the VGM. Finally, it replaces the
value with a mode picked according to the computed probabilities and a scalar
representing the value within the mode.

A CTGAN also integrates a conditional generator into the architecture of a GAN
to manage the class imbalances of categorical columns . It does so by
leveraging three key elements: the conditional vector, the generator loss, and
the training-by-sampling method. The conditional vector indicates if a
particular discrete column $D_{i^{*}}$ must have a particular value $k^{*}$.
Since the $D_{i^{*}}$ is represented as a one-hot encoding $d_{i^{*}}$, the
condition can be written $d_{i^{*}}^{(k^{*})}=1$. Each discrete column $D_{i}$
has a corresponding mask vector $m_{i}$ shaped like its one-hot encoding
$d_{i}$ to represent the condition

$m_{i}^{(k)}=\begin{cases}1&\text{if }i=i^{*}\text{ and }k=k^{*}\\\
0&\text{otherwise}.\end{cases}$

The conditional vector is the concatenation of these masks. During the feed-
forward pass, nothing encourages the conditional generator to produce the
desired output  $D_{i^{*}}=k^{*}$.  To enforce it, the conditional generator
loss is penalized by adding the cross-entropy between the generated discrete
columns and the conditional vector. The training-by-sampling approach samples
the conditional vector and training data to help the model evenly explore all
possible values in discrete columns. It first randomly selects a discrete
column, then constructs a probability mass function across the range of values
for that column, and samples a value according to this mass function. Finally,
it creates the conditional vector representing the sampled column and value.
We refer to Xu et al. [60] for a deeper dive into CTGAN’s formalism.

###  3.3 Tabular variational autoencoder

The architecture of a variational autoencoder (VAE) [37], considered by
Borysov et al. [13] in the context of synthetic populations, is analogous to
an autoencoder. A probabilistic encoder $q_{\phi}(Z\mid X)$ is meant to map
the input to a multivariate latent distribution and a probabilistic decoder
$p_{\theta}(X\mid Z)$ maps back the latent distribution to the data space. The
probabilistic encoder and decoder are represented by neural networks, with
parameters $\phi$ and $\theta$, respectively.

We would like to approximate the data distribution by first sampling from a
known prior distribution. Let $z$ represent a latent encoding of a data point
$x$ and consider their joint distribution $p_{\theta}(X,Z)$. Conditioning on
$Z$ gives us the parameterized approximation $p_{\theta}(X)$ of the data
distribution:

$p_{\theta}(x)=\int_{z}p_{\theta}(x\mid z)p_{\theta}(z)dz.$

The computation of $p_{\theta}(x)$ is usually intractable, for instance  when
likelihood functions $p_{\theta}(x\mid z)$ correspond to neural networks with
nonlinear hidden layers. For this reason,  it is necessary to introduce an
auxiliary parameterized family of functions, typically the  Normal
distributions ${\mathcal{N}}(\mu(x),\Sigma(x))$, to approximate the posterior
distribution

$q_{\phi}(z\mid x)\approx p_{\theta}(z\mid x).$

as schematized in Figure 2.  The objective is to make $q_{\phi}(z\mid x)$ as
close as possible to $p_{\theta}(z\mid x)$ by minimizing a distance measure
between distributions. In addition, we would like to optimize the generative
model to reduce the reconstruction error between the input and output.
Maximizing the evidence lower bound (ELBO) defined as

$\displaystyle L(x)$ $\displaystyle={\mathbb{E}}_{z\sim q_{\phi}(\cdot\mid
x)}\left[\frac{p_{\theta}(x,z)}{q(z\mid x)}\right]$ $\displaystyle=\ln
p_{\theta}(x)-D_{KL}\left(q(z\mid x)\parallel p_{\theta}(z\mid x)\right)$

is equivalent to minimizing the reverse Kullback-Leibler divergence between
the approximation and true posteriors and jointly maximizing the log-
likelihood of the data.  This is typically achieved by using a stochastic
gradient descent procedure. For more details,  we refer the reader to [37].

Figure 2: Schematic description of a variational autoencoder for which
$q_{\phi}(\cdot\mid x)$ is assumed to be the family of  Normal distributions.

A tabular variational autoencoder (TVAE) [60] adapts a VAE to tabular data
with the same preprocessing and loss function modification procedures as in
the CTGAN framework.

###  3.4 Copula-based generation

Most of household surveys and similar studies rely on detailed but small
samples, while we can access to more complete information , not captured in
such surveys, regarding specific  household characteristics. In particular,
data from household surveys are typically presented as multivariate vectors
implicitly capturing the dependencies between their components, but we have
sometimes access to the distributions of individual marginals at the
population level. Another typical situation is when we want to focus on a sub-
area covered by the sample, resulting in a limited number of observation
vectors, while we have a comprehensive knowledge of the distribution of the
individual factors. Finally, we sometimes  work with areas where we have
access to the marginal distributions only, and no household vector data. For
example, in a census data generation context, it is reasonable to conjecture
that marginals encode demographic elements and that multidimensional
dependencies encode more complex socioeconomic patterns shared across related
regions. Assuming that the dependencies can be adequately captured by a sample
over a different geographical sector, we aim to be able to exploit this
dependencies structure over the area of interest. We thus need a tool to
separate the dependencies structure to the specific marginal distributions,
which is provided by the copula theory.

A $d$-dimensional copula $C$ is a multivariate cumulative distribution
function (CDF) on $[0,1]^{d}$ having all marginals uniformly distributed on
$[0,1]$ [43, 33, 34, 45]. A fundamental result, due to Sklar [56], states that
any multivariate distribution can be represented by means of a copula.

######  Theorem (Sklar’s theorem).

Let $H$ be a multivariate distribution function with margins
$F_{1},\ldots,F_{d}$, then there exists a copula $C$ such that

$H(x_{1},\ldots,x_{d})=C\left(F_{1}(x_{1}),\ldots,F_{d}(x_{d})\right),\quad
x_{1},\ldots,x_{d}\in\overline{{\mathbb{R}}}.$ (1)

If $F_{j}$, $j=1,\ldots,d$, are continuous then $C$ is unique. Otherwise $C$
is uniquely determined on the Cartesian product of the range of the marginals
$F_{1}\left(\overline{{\mathbb{R}}}\right)\times\ldots\times
F_{d}\left(\overline{{\mathbb{R}}}\right)$. Conversely, if $C$ is a copula and
$F_{1},\ldots,F_{d}$ are univariate distribution functions, then  the function
$H$ defined above is a multivariate distribution function with margins
$F_{1},\ldots,F_{d}$.

Let $X$ be the random vector behind the population ${\mathcal{X}}$, from which
we extract a sample, but we are interested in the population ${\mathcal{Y}}$,
governed by the random vector $Y$. We assume from now that $X$ and $Y$ share
the copula $C$. This assumption is obviously valid when ${\mathcal{X}}$ is a
Monte Carlo sample from ${\mathcal{Y}}$, and we argue in the following that it
also holds when one population is a subset from the second one or when they
are two different populations with similar characteristics. We formalize this
idea by introducing the definition of a shared copula.

######  Definition (Shared copula).

The random vectors $X$ and $Y$ are said to share a copula $C$ if there exists
a copula $C$ such that for any vector
$(x_{1},\ldots,x_{d})\in\overline{{\mathbb{R}}}^{d}$,
$F^{X}(x_{1},\ldots,x_{d})=C\left(F^{X}_{1}(x_{1}),\ldots,F^{X}_{d}(x_{d})\right)$
and
$F^{Y}(x_{1},\ldots,x_{d})=C\left(F^{Y}_{1}(x_{1}),\ldots,F^{Y}_{d}(x_{d})\right)$.

Since a copula is a multivariate distribution on $[0,1]^{d}$, the first step
is to cast the observations as vectors in $[0,1]^{d}$. This can be easily
achieved as it is standard in machine learning to apply some normalization
procedure to the data (see for instance Larose and Larose [40]), for instance
by considering the empirical CDF (ECDF) of each of the feature $X_{i}$,
$i=1,\ldots,d$, defined as

$\hat{F}_{i}(x)=\frac{1}{N}\sum_{i=1}^{N}\mathds{1}_{x_{i}\leq x},$ (2)

where $\mathds{1}$ is the indicator function. The multivariate distribution
underlying the normalized observations can be estimated using any of the
previous ly considered methods.

Suppose now we wish to generate synthetic data of a target population from
which we only have marginal information and that we have access to a sample of
another, source, population sharing the structure of the target. Our
population synthesis procedure is summarized in Algorithm 1.

  1. Step 1

Normalize the source population data using the ECDFs $\hat{F}_{i}(\cdot)$,
$i=1,\ldots,d$.

  2. Step 2

Train the model on the normalized data to learn a copula $C$.

  3. Step 3

Generate a synthetic population of vectors in $[0,1]^{d}$ by sampling from
$C$.

  4. Step 4

Transform any generated vector ${\boldsymbol{u}}=(u_{1},\ldots,u_{d})$ in a
vector ${\boldsymbol{y}}$ in the target population as

${\boldsymbol{y}}=\left((F_{1}^{Y})^{-1}(u_{1}),\ldots,(F_{d}^{Y})^{-1}(u_{d})\right),$

where $(F^{Y}_{i})^{-1}(\cdot)$ is the pseudo-inverse distribution function of
the $i$-th target marginal, defined as

$(F^{Y}_{i})^{-1}(u)=\min\left\\{x_{i}:F^{Y}_{i}(x_{i})\geq u\right\\}.$

Algorithm 1 Synthetic population generation

Algorithm 1 implicitly assumes that the marginal distribution functions of the
source and target populations have the same ranges, an assumption often
violated, especially when the marginals follow discrete distributions, which
is especially the case when we use the ECDFs (2), as illustrated in Figure
3(a). We address this issue by relaxing any discrete distribution function as
a continuous distribution with the following heuristic. Consider the ECDF
$\hat{F}_{i}(\cdot)$, with the range
$\\{x_{i}^{(1)},\ldots,x_{i}^{(m_{i})}\\}$, and $x_{i}^{(p)}<x_{i}^{(q)}$ if
$p<q$. We construct the relaxed ECDF $\tilde{F}_{i}(\cdot)$ as the continuous
piecewise linear function obtained by considering the linear interpolation s
between consecutive values
$\hat{F_{i}}\left(x_{i}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(k)}}\right)$
and
$\hat{F_{i}}\left(x_{i}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(k+1)}}\right)$,
$k=1,2,\ldots,m_{i}-1$, as depicted in Figure 3(b). The relaxed marginals
being now continuous, each one is uniformly distributed on $[0,1]$. We now
extend the copula $C$ to the domain $[0,1]^{d}$ by setting
$C\left(\tilde{F}_{1}(x_{1}),\ldots,\tilde{F}_{d}(x_{d})\right)$ as the linear
interpolation, component by component,  from
$C\left(\hat{F}_{1}\left(x_{1}^{(k_{1})}\right),\ldots,\hat{F}_{i}\left(x_{i}^{(k_{i})}\right),\ldots,\hat{F}_{d}\left(x_{d}^{(k_{d})}\right)\right)$
to
$C\left(\hat{F}_{1}\left(x_{1}^{(k_{1})}\right),\ldots,\hat{F}_{i}\left(x_{i}^{(k_{i}+1)}\right),\ldots,\hat{F}_{d}\left(x_{d}^{(k_{d})}\right)\right)$,
with $x_{i}^{(k_{i})}\leq x_{i}\leq x_{i}^{(k_{i}+1)}$,
$k=1,2,\ldots,m_{i}-1$, $i=1,\ldots,d$. Note that the copula built on the
relaxed distribution functions still satisfies Sklar’s theorem for the source
marginals since it produces the same realizations on the Cartesian product of
the range of the source marginals, but is now uniquely defined as any marginal
distribution $\tilde{F}_{i}(\cdot)$, $i=1,\ldots,d$ follows an uniform
distribution on $[0,1]$.

$0$$0.2$$0.4$$0.6$$0.8$$1$$x_{i}^{(1)}$$x_{i}^{(2)}$$x_{i}^{(3)}$$x_{i}^{(4)}$$\hat{F}_{i}(x)$
(a) Empirical cumulative distribution function.

$0$$0.2$$0.4$$0.6$$0.8$$1$$x_{i}^{(1)}$$x_{i}^{(2)}$$x_{i}^{(3)}$$x_{i}^{(4)}$$\tilde{F}_{i}(x)$
(b) Empirical cumulative distribution function with piecewise linear
interpolations.

Figure 3:

In short, we replace Step 3 in Algorithm 1 by

  1. Step 3

Generate a synthetic population of vectors in $[0,1]^{d}$ by sampling from $C$
extended to $\tilde{F}_{i}(\cdot)$, $i=1,\ldots,d$.

##  4 Numerical experiments

To evaluate the methods described in the previous section, we generate
synthetic populations at the state, county, and census tract levels. For that
purpose, we use a sample of variables obtained from the American Community
Survey (ACS), as well as their actual distributions (marginals), collected
from the Decennial Census of Population Housing Data and the Internal Revenue
Service (IRS). The ACS is a national demographic survey conducted every year
by the U.S. Census Bureau that gathers sociodemographic and socioeconomic
information at the household and individual level. This information is
provided for geographical units called Public Use Micro Areas (PUMAs) that
contain at least a hundred thousand individuals and do not overlap or nest
within a single state. In our case, we use the 5-year ACS sample for the years
2012 to 2016 to ensure having a representative sample of the American
population, as well as data from the 2010 census for the marginals. Table 1
below lists the nine selected variables, definitions, level of aggregation,
and associated values. With respect to the marginals, the decennial census
data is updated every ten years at years ending with zeros and provides actual
counts of the total number of individuals and households for several
demographics and socioeconomic variables, including those of interest in our
study.

Table 1: List of variables. Name |  Definition |  Level |  Values
---|---|---|---
AGEP |  Age of person |  Individual |  0 to 99 years
SEX |  Gender of person |  Individual |  1: Male, 2: Female
RAC1P |  Race of person |  Individual |  1: White alone, 2: Black or African American alone, 3:: American Indian alone, 4: Alaska Native alone 5: American Indian and Alaska Native and no other races,6: Asian alone, 7: Native Hawaiian and Other Pacific Islander alone, 8: Some other race alone, 9: Two or more major race groups
ESR |  Employment status |  Individual |  1: Civilian employed, at work, 2: Civilian employed, with a job but not at work, 3: Unemployed, 4: Armed Forces, At Work, 5: Armed Forces, With a Job but Not at Work, 6: Not in Labor Force
HINCP |  Household income (past 12 months) |  Household |  1: $1 to less $25k, 2: $25k to less $50k, 3: $50k to less $75k, 4: $75k to less $100k, 5: $100k to less $200k, 6: $200k or more
HHT |  Household/family type |  Household |  1: Married couple household, 2: Male householder, no spouse present, 3: Female householder, no spouse present, 4: Male householder: Living alone, 5: Male householder: Not living alone, 6: Female householder: Living alone, 7:Female householder: Not living alone
NP |  Number of persons the household |  Household |  1,6: Number of persons in household 7: Household of 7 persons or more
WIF |  Workers in family during the past 12 months |  Household |  0: No workers, 1: 1 worker, 2: 2 workers, 3: 3 or more workers in family
HUPAC |  HH presence and age of children |  Household |  1: With children under 6 years only, 2: With children 6 to 17 years only, 3: With children under 6 years and 6 to 17 years, 4: No children

In order to assess the quality and goodness of fit of the synthetic data we
report two metrics: _the standardized root mean squared error_ (SRMSE) and the
count known as _sampling zeros_ (SZ). Regarding the SRMSE, we compute it as in
Sun and Erath [58]:

$SRMSE=\sqrt{M\sum_{m_{1}=1}^{M_{1}}\ldots\sum_{m_{d}=1}^{M_{d}}\left(\pi_{m_{1}\ldots
m_{d}}-\hat{\pi}_{m_{1}\ldots m_{d}}\right)^{2}},$ (3)

where $\pi_{m_{1}\ldots m_{d}}$ and $\hat{\pi}_{m_{1}\ldots m_{d}}$ are the
relative frequencies of a particular combination in the reference data and in
the synthetic data, respectively and $M=\prod_{i=1}^{d}M_{i}$. SRMSE captures
whether a combination of synthetic data appears in the actual data in similar
proportions. A value of $0$ means a perfect match while larger values evidence
increasing distance between true and synthetic data. Moreover, we report the
average SRMSE over subsets of $n$ dimensions; when $n=1$, the SRMSE measures
the fitting of the marginals, while when $n>1$ it measures the fitting of the
multivariate dependencies.

SRMSE does not take into account the diversity of data. A synthetic
combination might be desired even if it does not appear in the reference data.
To assess the diversity of the produced synthetic data, we implement SZ, the
count of the combinations of variables that are in the test set but not in the
training set [23]. A SZ count of $0$ means that the generative model was
unable to produce unseen realistic examples that do not appear in the training
set, but do appear in the test set. Larger sampling zero values means that the
generative model can output out-of-sample examples.

In our experiments, we use the BN implementation from the Pomegranate package
for Python [54] and the CTGAN and TVAE implementations from The Synthetic Data
Vault [47].

###  4.1 Population synthesis for the State of Maryland

This section is devoted to the task of generating synthetic data for a target
population of which only a sub-sample and the marginal distributions are
available. We compare BN, CTGAN, and TVAE to their copula counterparts, which
means that we use them in the second and third steps of Algorithm 1 with the
objective of determining whether the copula framework improves their
performances in matching the marginal distributions of the target. These
methods for which we applied the copula framework are respectively referred to
as BN Copula, CTGAN Copula, and TVAE Copula.  For BN we use the heuristic
described in Section 3.1, while for CTGAN and TVAE we implement the network
structures and hyperparameters suggested by Xu et al. [60]. We also report the
performances of IPF and that of an independent baseline created from the
source population by sampling its variables independently.

We generate synthetic populations at the state and county levels. At the
geographical level considered, the target population is the aggregation of
available PUMAs, thus making it possible to assess the fitting of multivariate
dependencies. To mimic how public agencies generally provide a subset from the
whole population, we created a source population by randomly sampling $1\%$
from the target. In addition to the source population,  the empirical marginal
distributions of all the target population’s variables were also provided to
IPF, BN Copula, CTGAN Copula, and TVAE Copula.

As indicated above, in order to evaluate the accuracy of the relationships
between the variables of the generated data, we compute the average SRMSE
between the synthetic and target populations over the possible combinations of
$n$ variables for $n$ ranging from $1$ to the total number of dimensions. We
also report the number of sampled zeros to determine which models can generate
unseen realistic examples.

Table 2 presents SRMSE and SZ when the target population is the Sate of
Maryland. It can be observed that IPF best captures marginal, bivariate, and
trivariate distributions, but fails to capture higher order interactions
between variables. In accordance with the findings of Sun and Erath [58],
unlike IPF, BN is able to capture complex relationships between variables.
However, it is no better than the independent baseline at matching the
target’s marginals, which motivates the use of the copula framework. BN Copula
indeed improves the SRMSE score for marginal distributions compared to BN.
Overall, IPF performs best in small dimensions, closely followed by BN Copula,
but, as dimensions increase, BN Copula significantly outperforms IPF.

| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | SZ
---|---|---|---|---|---|---|---|---|---|---|---
Ind | 0.049 | 0.433 | 1.174 | 2.462 | 4.709 | 8.732 | 16.166 | 30.503 | 59.633 | 122.089 | 6937
CTGAN | 0.399 | 0.866 | 1.630 | 2.976 | 5.438 | 10.067 | 19.007 | 36.761 | 73.087 | 149.743 | 23053
CTGAN Copula | 0.362 | 0.789 | 1.490 | 2.723 | 4.965 | 9.154 | 17.199 | 33.136 | 65.856 | 135.591 | 24786
TVAE | 0.425 | 1.099 | 2.249 | 4.271 | 7.899 | 14.515 | 26.770 | 49.906 | 94.666 | 183.920 | 18792
TVAE Copula | 0.409 | 1.050 | 2.138 | 4.033 | 7.393 | 13.448 | 24.542 | 45.313 | 85.366 | 165.563 | 21230
BN | 0.047 | 0.222 | 0.591 | 1.326 | 2.779 | 5.680 | 11.586 | 23.983 | 50.957 | 111.703 | 42447
BN Copula | 0.009 | 0.183 | 0.546 | 1.270 | 2.710 | 5.599 | 11.494 | 23.881 | 50.841 | 111.555 | 41575
IPF | 0.001 | 0.129 | 0.496 | 1.567 | 4.592 | 12.925 | 35.395 | 94.869 | 249.543 | 644.951 | 0
Table 2: State-level SRMSE values between the target population and synthetic
agents . Column $n$ represents the SRMSE averaged over the possible $n$-tuples
variables. The last column lists sampling zero counts. Figure 4: Marginal fit
of the county for the state-level experiment.

The county-level experiment, whose target is the aggregation of PUMAs in Anne
Arundel County, shows very similar results. Table 3 presents the SRMSE and SZ
metrics; BN Copula consistently succeeds in capturing the marginal
distributions as well as the multivariate dependencies, outperforming IPF even
in small dimensions.

Figures 4 and 5 plot the marginal distributions of the County and PUMA
variables respectively from the  state and  county experiments. It can be
observed that, compared to BN, BN Copula’s counts are closer to those of the
target.

| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | SZ
---|---|---|---|---|---|---|---|---|---|---|---
Ind | 0.165 | 0.560 | 1.309 | 2.655 | 5.087 | 9.623 | 18.477 | 36.655 | 75.779 | 163.087 | 291
CTGAN | 0.361 | 0.865 | 1.688 | 3.103 | 5.622 | 10.270 | 19.257 | 37.588 | 76.894 | 164.451 | 459
CTGAN Copula | 0.268 | 0.688 | 1.421 | 2.725 | 5.102 | 9.578 | 18.368 | 36.478 | 75.530 | 162.715 | 480
TVAE | 0.271 | 0.770 | 1.786 | 3.848 | 8.035 | 16.554 | 33.941 | 69.598 | 143.166 | 295.982 | 923
TVAE Copula | 0.167 | 0.533 | 1.272 | 2.729 | 5.605 | 11.302 | 22.732 | 46.102 | 94.990 | 199.617 | 1216
BN | 0.165 | 0.506 | 1.163 | 2.388 | 4.670 | 9.024 | 17.673 | 35.665 | 74.711 | 162.174 | 672
BN Copula | 0.015 | 0.303 | 0.899 | 2.042 | 4.213 | 8.431 | 16.928 | 34.771 | 73.691 | 161.081 | 748
IPF | 0.031 | 0.540 | 1.972 | 5.827 | 15.804 | 41.012 | 103.596 | 256.709 | 626.142 | 1505.718 | 0
Table 3: County -level SRMSE values between the target population and
synthetic agents, column $n$ represents the SRMSE averaged over the possible
$n$-tuples variables. The last column lists sampling zero counts. Figure 5:
Marginal fit of the PUMA for the county-level experiment.

In both experiments, the expected inability of IPF to generalize to new
realistic out-of-sample examples is exposed. This is of great relevance since,
although it may be a good procedure to capture the marginal distributions and
the low-dimensional relationships, the synthetic populations it generates
suffer from a significant lack of diversity.

It can also be noticed that all methods in conjunction with copulas perform
better than their standalone counterparts. Moreover, a common trend among
these methods is that the  observed performance gaps narrow as the number of
dimensions increases. We argue that this advantage of the copula framework is
due to its ability to incorporate and match marginal distributions of the
target, and that the effect of this advantage is less and less present as
inter-dimensional relationships become more complex in size.

We should also mention that the proposed method heavily relies on the
underlying model: its performance will follow the capacity of the generative
model to learn the copula. CTGAN and TVAE poorly learn the population’s
distribution, which is why we can’t expect their copula counterparts to
perform well. Conversely, the proposed framework inherits the complexity and
diversity of the underlying model when it is well designed for the given
modeling task. According to these two state and county-level experiments, BN
Copula is the most consistent model; like classical BN, it captures high-
dimensional relationships between variables and is able to sample zero cells,
also matching the marginal distributions of the target.

###  4.2 Population synthesis at the census tract level

In this section we show the results of testing our methodology at a smaller
geographical level. We analyze the case when the source sample is a PUMA and
the target population is a census tract inside that PUMA. In this experiment
no sample of the target population is known. Therefore, classical ML
frameworks are not applicable. However, we do know the marginal distributions
of NP, RAC1P, SEX, HUPAC, and HHT at the census tract level from the decennial
census 2010. We also estimate the marginal distribution for HINCP at the
census tract level from the IRS data. We do so by counting the total number of
tax returns by the different income categories, which included the totals of
single returns, joint returns, head of family returns, and the number of
dependents. We use these marginal distributions to estimate a BN Copula at
PUMA level and generate the synthetic population at census tract level
assuming that the dependency structure can be maintained across the two
geographical levels. We compare these results with those obtained with the
modified IPF implemented in PopGen  2.0 [41]; this procedure uses two level of
spatial resolutions and, in particular, applies constraints at both PUMA and
census tract levels. SRMSE cannot be computed in this case, and the two
methods are compared based on the marginals for the variables considered.
Results relative to PUMA 1201 and census tract 702204 in Anne Arundel County
are shown in Figure 6.

At this small geographical level, IPF is close to the source but fails to
match the marginals of the target for most of the variables included in our
analysis. We argue that IPF might  be fac ing the zero-cell problem and ,
therefore, it is unable to  recreate the observations that are in the target
population [28] . Namely, when a demographic group of the target population is
not represented in the source sample, its weights remain zero in the
contingency table so that it cannot be sampled at the allocation stage. In
contrast, the BN Copula method is able to produce a significant number of out-
of-sample observations and fits almost perfectly the target’s marginals. These
results allow us to conclude that the proposed copula framework combined with
BN is able to reproduce marginals at small-area level. This is relevant for
producing statistics and estimates at small -area level, which are of
increasing interest to researchers and planning agencies.

Figure 6: Marginals fit for tract 702204 when source is PUMA 1201

###  4.3 Learning the dependencies from one PUMA to another

In previous experiments at the state, county, and census tract levels, the
source population was obtained from a subregion of the target population.
However, as stated in Section 3, we claim that different populations with
similar characteristics may share a copula even if the source is not  a subset
of the target. In this experiment, we propose a spatial transferability , from
one PUMA to another, of a trained model that captures the dependencies of the
selected variables. The concept of model transferability has been widely used
in transportation forecasting modeling between different regions when  no or
partial data is available, which reduces the time, cost, and effort of the
data collection process [55, 14].  To some extent, Goulias et al. [26]
performed a similar study, using the IPF procedure to generate a synthetic
population carrying information from county to census tract level, to estimate
the total number of daily trips per person.

We  tested our idea on several pairs of PUMAs and they all give similar
results. Here we exhibit the results for a random pair of PUMAs selected to be
the source and target populations.

The model transfer procedure cannot be executed on the classical ML methods,
therefore we compare IPF and the copula-based frameworks since they both can
incorporate the knowledge of the target’s marginal distributions. According to
the findings of our experiments at the state and county levels, BN seems more
suited to the relatively low-dimensional synthetic population generation
problem at hand, hence we report in Table 4 the performance of the copula
framework only applied to the BN. The SRMSE values show that that IPF is
slightly better than BN Copula at fitting the marginals for low-dimensional
settings (up to five dimensions). However, as the dependencies grow, IPF again
fails to capture the relationships while the SRMSE for the copula framework
stays below the independent baseline. It should be stressed that IPF only
replicates observations from the source and will always show null values in
the sampling zeros column.

| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | SZ
---|---|---|---|---|---|---|---|---|---|---
Ind | 0.179 | 0.634 | 1.601 | 3.540 | 7.513 | 15.971 | 34.577 | 76.500 | 172.372 | 37
BN Copula | 0.040 | 0.281 | 0.901 | 2.342 | 5.654 | 13.289 | 30.958 | 71.979 | 167.354 | 182
IPF | 0.042 | 0.218 | 0.722 | 2.033 | 5.341 | 13.577 | 33.818 | 82.975 | 200.853 | 0
Table 4: SRMSE values between the target population and synthetic agents for
the PUMA to PUMA experiment . Column $n$ represents the SRMSE averaged over
the possible $n$-tuples variables. The last column lists sampling zero counts.

##  5 Conclusion

In this paper, we introduce a framework that integrates the theory of copulas
into traditional machine learning generative methods to capture both joint and
marginal distributions in synthetic population generation. We also explore the
possibility to transfer the copula trained on one region and learn about the
joint distribution of the variables in a different one. The classical data
normalization relying on empirical distribution functions of the features
allows us to consider the observations as realizations of an unknown copula,
capturing the dependencies between their components, which can be learned
using standard techniques. We propose a heuristic approach to overcome the
issues related to discrete distributions, and propose a simple algorithm to
generate a synthetic population with different marginals but the same
dependencies structure. Since the geographical level is a crucial factor for
population synthesis and can influence the performance of different deep
generative methods, we conduct comparison experiments at the state, county,
and census tract levels. We apply the proposed copula framework to each of the
selected ML generative methods: BN, CTGAN, and TVAE, which are suitable for
the size of our problem and the nature of our data. Both accuracy and
diversity are important criteria to evaluate the performance of different
synthetic population generation methods, so we use SRMSE to evaluate the
accuracy of the synthetic population and sampling zeros to evaluate the
diversity. We compare the performance of these methods and their copula
counterparts to IPF and a baseline population obtained from the ACS sample
source by bootstrapping its dimensions independently.

The results of the experiments at the three geographical levels considered
show that the copula framework improves the performance of all the ML methods
under study. Also, the BN Coppula always outperforms the BN method at both
state and county level s based on the SRMSE metric and the ability to recover
the marginal distributions of the variables that characterize the population.
At census tract level, where no sample is available , making conventional ML
techniques not applicable, the IPF with double geographical constraints fails
to recover the target values for most of the variables considered, while the
BN Copula follows the expected patterns. The proposed spatial transferability,
which is again not feasible with standard ML techniques,  not only shows that
at PUMA level the BN Copula framework always produces results that are
superior to the independent case, but that they are close, or are even better,
to  those delivered by IPF (up to five dimensions, or more than five,
respectively) , while  also generating observations that are not in the
original sample.

Nevertheless, it is worth mentioning that our study included only discrete
variables and not continuous  ones. Hence, we have not been able to explore
the performance of the proposed copula framework in these cases. This is an
extension of this work that we plan to carry out in the future. Namely, to
explore more ML generative methods and different data types to prove our
conclusions in different methodological and information contexts. Finally,
this study opens up to the opportunity to combine data from multiple sources
(e.g. disaggregated survey data and aggregated administrative data) to
overcome the limited sample size problem for small areas. This is especially
important when producing population statistics and travel indicator estimates
at small geographical level, which are critical for policy analysis in
governmental agencies.

## Acknowledgments

The work of Fabian Bastin is supported by the Natural Sciences and Engineering
Research Council of Canada [Discovery Grant 2022-04400]. Pascal Jutras-Dubé
was furthermore supported by a graduate grant under the NSERC CREATE Program
on Machine Learning in Quantitative Finance and Business Analytics (Fin-ML).

## References

  * Arentze and Timmermans [2004] T. A. Arentze and H. J. Timmermans.  A learning based transportation oriented simulation system.  _Transportation Research Part B_ , 38(7):613–633, 2004.
  * Arentze et al. [2007] T. A. Arentze, H. J. P. Timmermans, and F. Hofman.  Creating synthetic household populations: Problems and approach.  _Transportation Research Record_ , 2014:85–91, 2007.
  * Auld et al. [2012] J. Auld, M. Javanmardi, and A. Mohammadian.  Integration of activity scheduling and traffic assignment in ADAPTS activity-based model.  In _TRB 91st Annual Meeting Compendium of Papers DVD_ , number 12-4225, Washington DC, USA, Jan. 2012. Transportation Research Board.
  * Auld et al. [2009] J. A. Auld, A. Mohammadian, and K. Wies.  Population synthesis with subregion-level control variable aggregation.  _Journal of Transportation Engineering_ , 135(9):632–639, 2009.
  * Avramidis et al. [2009] A. Avramidis, N. Channouf, and P. L’Ecuyer.  Efficient correlation matching for fitting discrete multivariate distributions with arbitrary marginals and normal-copula dependence.  _INFORMS Journal on Computing_ , 21:88–106, 2009.
  * Barthelemy and Toint [2013] J. Barthelemy and P. L. Toint.  Synthetic population generation without a sample.  _Transportation Science_ , 47(2):266–279, 2013\.
  * Bhat and Eluru [2009] C. R. Bhat and N. Eluru.  A copula-based approach to accommodate residential self-selection effects in travel behavior modeling.  _Transportation Research Part B_ , 43(7):749–765, 2009.
  * Bishop [2006] C. M. Bishop.  _Pattern Recognition and Machine Learning_.  Springer, New-York, NY, USA, 2006.
  * Boquet et al. [2020] G. Boquet, A. Morell, J. Serrano, and J. L. Vicario.  A variational autoencoder solution for road traffic forecasting systems: Missing data imputation, dimension reduction, model selection and anomaly detection.  _Transportation Research Part C_ , 115:102622, 2020.
  * Borgomeo et al. [2015] E. Borgomeo, G. Pflug, J. W. Hall, and S. Hochrainer-Stigler.  Assessing water resource system vulnerability to unprecedented hydrological drought using copulas to characterize drought duration and deficit.  _Water Resources Research_ , 51(11):8927–8948, 2015.
  * Born et al. [2014] K. Born, S. Yasmin, D. You, N. Eluru, C. R. Bhat, and R. M. Pendyala.  Joint model of weekend discretionary activity participation and episode duration.  _Transportation Research Record_ , 2413:34–44, 2014.
  * Borysov and Rich [2021] S. S. Borysov and J. Rich.  Introducing synthetic pseudo panels: application to transport behaviour dynamics.  _Transportation_ , 48(5):2493–2520, 2021.
  * Borysov et al. [2019] S. S. Borysov, J. Rich, and F. C. Pereira.  How to generate micro-agents? A deep generative modeling approach to population synthesis.  _Transportation Research Part C_ , 106:73–97, 2019.
  * Bowman et al. [2014] J. L. Bowman, M. Bradley, J. Castiglione, and S. L. Yoder.  Making advanced travel forecasting models affordable through model transferability.  In _Transportation Research Board 93rd Annual Meeting_ , Jan. 2014\.
  * Bradley et al. [2010] M. Bradley, J. L. Bowman, and B. Griesenbeck.  Sacsim: An applied activity-based model system with fine-level spatial and temporal resolution.  _Journal of Choice Modelling_ , 3(1):5–31, 2010\.
  * Casati et al. [2015] D. Casati, K. Müller, P. J. Fourie, A. Erath, and K. W. Axhausen.  Synthetic population generation by combining a hierarchical, simulation-based approach with reweighting by generalized raking.  _Transportation Research Record_ , 2493:107–116, 2015.
  * Cherubini et al. [2004] U. Cherubini, E. Luciano, and W. Vecchiato.  _Copula Methods in Finance_.  John Wiley & Sons, Chichester, United Kingdom, 2004.
  * Cooper [1990] G. F. Cooper.  The computational complexity of probabilistic inference using Bayesian belief networks.  _Artificial Intelligence_ , 42(2):393–405, 1990\.
  * Deming and Stephan [1940] W. E. Deming and F. F. Stephan.  On a least squares adjustment of a sampled frequency table when the expected marginal totals are known.  _The Annals of Mathematical Statistics_ , 11(4):427–444, 1940.
  * Duguay et al. [1976] G. Duguay, W. Jung, and D. L. McFadden.  SYNSAM: A methodology for synthesizing household transportation survey data.  Working paper 7618, Institute of Transportation Studies, University of California, Berkeley, CA, USA, 1976.
  * Eluru et al. [2008] N. Eluru, A. R. Pinjari, J. Y. Guo, I. N. Sener, S. Srinivasan, R. B. Copperman, and C. R. Bhat.  Population updating system structures and models embedded in the comprehensive econometric microsimulator for urban systems.  _Transportation Research Record_ , 2076:171–182, 2008.
  * Farooq et al. [2013] B. Farooq, M. Bierlaire, R. Hurtubia, and G. Flötteröd.  Simulation based population synthesis.  _Transportation Research Part B_ , 58:243–263, 2013.
  * Garrido et al. [2020] S. Garrido, S. S. Borysov, F. C. Pereira, and J. Rich.  Prediction of rare feature combinations in population synthesis: Application of deep generative modelling.  _Transportation Research Part C_ , 120:102787, 2020.
  * Genest and Favre [2007] C. Genest and A.-C. Favre.  Everything you always wanted to know about copula modeling but were afraid to ask.  _Journal of hydrologic engineering_ , 12(4):347–368, 2007.
  * Goodfellow et al. [2014] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio.  Generative adversarial nets.  In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K. Weinberger, editors, _Advances in Neural Information Processing Systems_ , volume 27. Curran Associates, Inc., 2014.
  * Goulias et al. [2014] K. G. Goulias, S. K. Ravulaparthy, K. C. Konduri, and R. M. Pendyala.  Using synthetic population generation to replace sample and 3 expansion weights in household surveys for small area estimation of population parameters.  In _TRB 93rd Annual Meeting Compendium of Papers_ , number 14-0501, Washington DC, USA, Jan. 2014. Transportation Research Board.
  * Günthermann et al. [2020] L. Günthermann, I. Simpson, and D. Roggen.  Smartphone location identification and transport mode recognition using an ensemble of generative adversarial networks.  In _Adjunct Proceedings of the 2020 ACM International Joint Conference on Pervasive and Ubiquitous Computing and Proceedings of the 2020 ACM International Symposium on Wearable Computers_ , UbiComp-ISWC ’20, pages 311––316, New York, NY, USA, Sept. 2020. Association for Computing Machinery.
  * Guo and Bhat [2007] J. Y. Guo and C. R. Bhat.  Population synthesis for microsimulating travel behavior.  _Transportation Research Record_ , 2014:92–101, 2007.
  * Heckerman et al. [1995] D. Heckerman, D. Geiger, and D. M. Chickering.  Learning Bayesian networks: The combination of knowledge and statistical data.  _Machine Learning_ , 20(3):194–243, 1995.
  * Huang and Williamson [2001] Z. Huang and P. Williamson.  A comparison of synthetic reconstruction and combinatorial optimisation approaches to the creation of small-area microdata.  Working Paper 2001/2, University of Liverpool, Oct. 2001.
  * Hörl and Balac [2021] S. Hörl and M. Balac.  Synthetic population and travel demand for Paris and Île-de-France based on open and publicly available data.  _Transportation Research Part C_ , 130:103291, 2021.
  * Jaoua et al. [2013] A. Jaoua, P. L’Ecuyer, and L. Delorme.  Call-type dependence in multiskill call centers.  _Simulation_ , 89(6):722–734, 2013.
  * Joe [1997] H. Joe.  _Multivariate Models and Dependence Concepts_.  Springer, New York, NY, USA, 1997.
  * Joe [2015] H. Joe.  _Dependence modeling with copulas_.  CRC Press, Boca Raton, FL, USA, 2015.
  * Kao et al. [2012] S.-C. Kao, H. K. Kim, C. Liu, X. Cui, and B. L. Bhaduri.  Dependence-preserving approach to synthesizing household characteristics.  _Transportation Research Record_ , 2302:192–200, 2012.
  * Kaushik et al. [2019] K. Kaushik, C. Cirillo, and F. Bastin.  On modelling human population characteristics with copulas.  _Procedia Computer Science_ , 151:210–217, 2019.  The 10th International Conference on Ambient Systems, Networks and Technologies (ANT 2019) / The 2nd International Conference on Emerging Data and Industry 4.0 (EDI40 2019) / Affiliated Workshops.
  * Kingma and Welling [2014] D. P. Kingma and M. Welling.  Auto-encoding variational bayes.  In _International Conference on Learning Representations (ICLR)_ , Banff, AB, Canada, Apr. 2014.
  * Konduri et al. [2016] K. Konduri, D. You, V. Garikapati, and R. Pendyala.  Enhanced synthetic population generator that accommodates control variables at multiple geographic resolutions.  _Transportation Research Record_ , 2563:40–50, 2016.
  * Lam and Bacchus [1994] W. Lam and F. Bacchus.  Learning bayesian belief networks: An approach based on the mdl principle.  _Computational Intelligence_ , 10(3):269–293, 1994.
  * Larose and Larose [2014] D. T. Larose and C. D. Larose.  _Data Preprocessing_ , chapter 2, pages 16–50.  John Wiley & Sons, Hoboken, NJ, USA, 2014.
  * [41] Mobility Analytics Research Group.  MARG (2016) PopGen: Synthetic population generator.  http://www.mobilityanalytics.org/popgen.html.  Accessed December 4th, 2022.
  * Müller and Axhausen [2011] K. Müller and K. W. Axhausen.  Population synthesis for microsimulation: State of the art.  In _TRB 90th Annual Meeting Compendium of Papers DVD_ , number 11-1789, Washington DC, USA, Jan. 2011. Transportation Research Board.
  * Nelsen [2006] R. B. Nelsen.  _An Introduction to Copulas_.  Springer, New York, NY, USA, second edition, 2006.
  * Nikoloulopoulos [2013] A. K. Nikoloulopoulos.  Copula-based models for multivariate discrete response data.  In P. Jaworski, F. Durante, and W. K. Härdle, editors, _Copulae in Mathematical and Quantitative Finance_ , volume 213 of _Lecture Notes in Statistics_ , pages 231–249. Springer, 2013.
  * Okhrin et al. [2017] O. Okhrin, A. Ristig, and Y.-F. Xu.  _Copulae in High Dimensions: An Introduction_ , chapter 13, pages 247–277.  Springer, Berlin, Germany, third edition, 2017.
  * Oreshkin et al. [2016] B. Oreshkin, N. Réegnard, and P. L’Ecuyer.  Rate-based daily arrival process models with application to call centers.  _Operations Research_ , 64, 2016.
  * Patki et al. [2016] N. Patki, R. Wedge, and K. Veeramachaneni.  The synthetic data vault.  In _2016 IEEE International Conference on Data Science and Advanced Analytics (DSAA)_ , pages 399–410, Oct 2016.
  * Pinjari et al. [2009] A. R. Pinjari, C. R. Bhat, and D. A. Hensher.  Residential self-selection effects in an activity time-use behavior model.  _Transportation Research Part B_ , 43(7):729–748, 2009.
  * Pritchard and Miller [2012] D. R. Pritchard and E. J. Miller.  Advances in population synthesis: fitting many attributes per agent and fitting to household and person margins simultaneously.  _Transportation_ , 39(3):685–704, 2012.
  * Rana et al. [2010] T. A. Rana, S. Sikder, and A. R. Pinjari.  Copula-based method for addressing endogeneity in models of severity of traffic crash injuries: Application to two-vehicle crashes.  _Transportation Research Record_ , 2147:75–87, 2010.
  * Saadi et al. [2016a] I. Saadi, A. Mustafa, J. Teller, and M. Cools.  Forecasting travel behavior using Markov Chains-based approaches.  _Transportation Research Part C_ , 69:402–417, 2016a.
  * Saadi et al. [2016b] I. Saadi, A. Mustafa, J. Teller, B. Farooq, and M. Cools.  Hidden Markov model-based population synthesis.  _Transportation Research Part B_ , 90:1–21, 2016b.
  * Salvini and Miller [2005] P. Salvini and E. Miller.  ILUTE: An operational prototype of a comprehensive microsimulation model of urban systems.  _Networks and Spatial Economics_ , 5(2):217–234, 2005.
  * Schreiber [2018] J. Schreiber.  pomegranate: Fast and flexible probabilistic modeling in python.  _Journal of Machine Learning Research_ , 18(164):1–6, 2018.
  * Sikder et al. [2013] S. Sikder, A. R. Pinjari, S. Srinivasan, and N. Roosbeh.  Spatial transferability of travel forecasting models: a review and synthesis.  _International Journal of Advances in Engineering Sciences and Applied Mathematics_ , 5(2):104–128, 2013.
  * Sklar [1959] A. Sklar.  Fonctions de répartition à n dimensions et leurs marges.  _Publications de l’Institut de Statistique de l’Université de Paris_ , 8:229–231, 1959.
  * Smith et al. [1995] L. Smith, R. Beckman, K. Baggerly, D. Anson, and M. Williams.  TRANSIMS: Transportation analysis and simulation system.  techreport LA-UR-95-1641, Los Alamos National Laboratory, Los Alamos, NM, USA, July 1995.
  * Sun and Erath [2015] L. Sun and A. Erath.  A Bayesian network approach for population synthesis.  _Transportation Research Part C_ , 61:49–62, 2015.
  * Sun et al. [2018] L. Sun, A. Erath, and M. Cai.  A hierarchical mixture modeling framework for population synthesis.  _Transportation Research Part B_ , 114:199–212, 2018.
  * Xu et al. [2019] L. Xu, M. Skoularidou, A. Cuesta-Infante, and K. Veeramachaneni.  Modeling tabular data using conditional GAN.  In H. M. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché Buc, and E. B. Fox, editors, _Proceedings of the 33rd International Conference on Neural Information Processing Systems_ , pages 7335–7345, Red Hook, NY, USA, 2019. Curran Associates Inc.
  * Yazdizadeh et al. [2021] A. Yazdizadeh, Z. Patterson, and B. Farooq.  Semi-supervised GANs to infer travel modes in GPS trajectories.  _Journal of Big Data Analytics in Transportation_ , 3(3):201–211, 2021.
  * Ye et al. [2009] X. Ye, K. Konduri, R. Pendyala, B. Sana, and P. Waddell.  Methodology to match distributions of both household and person attributes in generation of synthetic populations.  In _TRB 88th Annual Meeting Compendium of Papers DVD_ , number 09-2096, Washington DC, USA, Jan. 2009. Transportation Research Board.
  * Yin et al. [2018] M. Yin, M. Sheehan, S. Feygin, J.-F. Paiement, and A. Pozdnoukhov.  A generative model of urban activities from cellular data.  _IEEE Transactions on Intelligent Transportation Systems_ , 19(6):1682–1696, 2018.
  * Zhang et al. [2019] D. Zhang, J. Cao, S. Feygin, D. Tang, Z.-J. Shen, and A. Pozdnoukhov.  Connected population synthesis for transportation simulation.  _Transportation Research Part C_ , 103:1–16, 2019.
  * Ziemke et al. [2018] D. Ziemke, J. W. Joubert, and K. Nagel.  Accessibility in a post-apartheid city: Comparison of two approaches for accessibility computations.  _Networks and Spatial Economics_ , 18(2):241–271, 2018.