Fit and compare distributions
Question. How close is a model's output (or a second dataset) to the observed data, and which metric should I use to say so?
You will get a gravity model fitted to observed flows between zones, its predicted flows, and that comparison made with every evaluation metric of the library: error measures, information measures, distribution distances, common-part-of-commuters measures and a spatio-temporal distance.
Prerequisites. Flows between places and Generate mobility with models.
Data. Brightkite check-ins (Stanford SNAP, 2008 to 2010): the first 300,000 rows (about 4 MB of the 60 MB file) with fastmob.data.load_brightkite, users with at least 100 check-ins, inside the contiguous United States. Zones are a synthetic 10 by 10 grid. For the occupancy distance each check-in is treated as a 30-minute stay (an assumption). The notebook stops if the download fails.
Observed flows between zones
We build the observed origin-destination flows between 100 grid zones from consecutive check-ins, exactly as in the population notebook.
import warnings
import fastmob
import geopandas as gpd
import numpy as np
import pandas as pd
from fastmob import models
from fastmob.data import load_brightkite
from fastmob.measures import evaluation
from shapely.geometry import box
from statsmodels.tools.sm_exceptions import SingularMatrixWarning
raw = load_brightkite(max_rows=300_000).to_pandas()
raw = raw[~(((raw["lat"] == 0) & (raw["lng"] == 0)) | (raw["lat"].abs() >= 90))]
counts = raw.groupby("uid").size()
us = raw[raw["uid"].isin(counts[counts >= 100].index)]
n, (lat0, lat1, lng0, lng1) = 10, (25, 50, -125, -65)
us = (
us[us["lat"].between(lat0, lat1) & us["lng"].between(lng0, lng1)]
.sort_values(["uid", "datetime"])
.reset_index(drop=True)
)
grid = gpd.GeoDataFrame(
{"tile_id": [str(k) for k in range(n * n)]},
geometry=[
box(
lng0 + j * (lng1 - lng0) / n,
lat0 + i * (lat1 - lat0) / n,
lng0 + (j + 1) * (lng1 - lng0) / n,
lat0 + (i + 1) * (lat1 - lat0) / n,
)
for i in range(n)
for j in range(n)
],
crs="EPSG:4326",
)
observed = (
fastmob.TrajDataFrame(us).to_flowdataframe(grid, self_loops=False).to_pandas()
)
visits_per_tile = (
fastmob.TrajDataFrame(us)
.mapping(grid, remove_na=True)
.to_pandas()
.groupby("tile_id")
.size()
)
print(
len(us),
"check-ins;",
len(observed),
"zone pairs with flow;",
int(observed["flow"].sum()),
"trips",
)
Fit and run a gravity model
A gravity model predicts the flow between two zones from their relevance (here the number of check-ins in each zone) and the distance between them. We fit it to the observed flows (Gravity.fit estimates the distance exponent), then generate predicted flows. The model is singly constrained: it keeps each zone's total outflow equal to the observed one.
outflow = observed.groupby("origin")["flow"].sum()
zones = grid.assign(
relevance=grid["tile_id"].map(visits_per_tile).fillna(0).astype(float),
tot_outflow=grid["tile_id"].map(outflow).fillna(0).astype(float),
)
catalogue = fastmob.Locations.from_tessellation(zones)
model = models.Gravity(gravity_type="singly constrained")
with warnings.catch_warnings():
# The origin indicator columns are collinear by construction.
warnings.simplefilter("ignore", SingularMatrixWarning)
model.fit(
fastmob.FlowDataFrame(observed, locations=catalogue),
relevance_column="relevance",
)
predicted = model.generate(catalogue, out_format="flows").to_pandas()
print(
"distance exponent:",
np.round(model.deterrence_func_args, 2),
"|",
len(predicted),
"predicted zone pairs; total predicted",
round(predicted["flow"].sum()),
"vs observed",
int(observed["flow"].sum()),
)
The fitted distance exponent is about -1.78: predicted flow falls roughly with distance to the power 1.8. The singly constrained model keeps each zone's total outflow, so the predicted total equals the observed 13,835 trips, but spreads it over 5,402 zone pairs against 1,117 pairs with observed flow, because a gravity model gives every pair some flow. As in the previous notebook, statsmodels warns that the design matrix is rank-deficient (kept in the output above). The warning is structural: the intercept and the one-indicator-per-origin columns of the singly constrained fit are collinear, so it appears on any data; the distance exponent is still estimated.
Error metrics on matched flows
To compare predicted with observed flows we put them on the same list of zone pairs (a pair with no flow counts as 0). Then the error metrics: mse (mean squared error), rmse (its root, in trips), nrmse (rmse divided by the total number of observed trips), r_squared (share of the variance of the observed flows that the prediction explains) and max_error (the largest signed error, observed minus predicted).
pairs = observed.merge(
predicted, on=["origin", "destination"], how="outer", suffixes=("_obs", "_pred")
).fillna(0)
true, pred = pairs["flow_obs"].to_numpy(), pairs["flow_pred"].to_numpy()
errors = {
name: getattr(evaluation, name)(true, pred)
for name in ("mse", "rmse", "nrmse", "r_squared", "max_error")
}
print(len(pairs), "zone pairs")
pd.Series(errors).round(4)
5402 zone pairs
mse 134.5349
rmse 11.5989
nrmse 0.0008
r_squared 0.9238
max_error 488.7042
dtype: float64
R² is 0.92, so the prediction explains 92% of the variance of the observed flows over the 5,402 zone pairs; the root mean squared error is 11.6 trips per pair; the largest error is 489 trips (observed above predicted). The normalised error (0.0008) looks tiny only because it divides by the total of 13,835 trips: do not compare it with values from other datasets.
Information and distribution metrics
These compare shapes. After normalising both flow tables to sum to 1 we compute the Kullback-Leibler divergence and the information gain (two related measures of the information lost when the prediction replaces the observation, in nats; they are defined only where both are positive, so we use the zone pairs with observed flow and add a tiny constant to the prediction), the Jensen-Shannon divergence (symmetric, between 0 and ln 2 in nats) and the Wasserstein distance between the two sets of flow values. For all four, 0 means identical.
p_true, p_pred = true / true.sum(), pred / pred.sum()
eps = 1e-12
shapes = {
"kullback_leibler_divergence": evaluation.kullback_leibler_divergence(
p_true[true > 0], p_pred[true > 0] + eps
),
"information_gain": evaluation.information_gain(
p_true[true > 0], p_pred[true > 0] + eps
),
"jensen_shannon_divergence": evaluation.jensen_shannon_divergence(p_true, p_pred),
"wasserstein_distance": evaluation.wasserstein_distance(true, pred),
}
pd.Series(shapes).round(4)
kullback_leibler_divergence 0.3426
information_gain 0.4425
jensen_shannon_divergence 0.1069
wasserstein_distance 0.5476
dtype: float64
The Jensen-Shannon divergence is 0.107 (its maximum is 0.69), the Kullback-Leibler divergence 0.34, the information gain 0.44 and the Wasserstein distance 0.55 trips. These are moderate: the model reproduces the main pattern but not the detail. Compare each number only with the same metric for another model, not across metrics.
Common part of commuters and links
common_part_of_commuters (CPC) is the share of trips that observed and predicted flows send along the same links (1 identical, 0 disjoint); common_part_of_links (CPL) the share of links they have in common. common_part_of_commuters_distance (CPCD) compares the distance distribution of two sets of trips in 2 km bins, so it needs Trips with a distance_km column. For CPCD we use trips between H3 cells (resolution 7, centre to centre distance) of two random halves of the users; it is 1 for identical distance distributions.
def with_distance(trips):
table = trips.to_pandas()
cells = pd.unique(
pd.concat([table["origin_location_id"], table["destination_location_id"]])
)
centres = fastmob.h3_to_latlng(pd.DataFrame({"h3_cell": cells})).set_index(
"h3_cell"
)
a, b = (
centres.loc[table["origin_location_id"]].to_numpy(),
centres.loc[table["destination_location_id"]].to_numpy(),
)
phi1, phi2, dlng = (
np.radians(a[:, 0]),
np.radians(b[:, 0]),
np.radians(b[:, 1] - a[:, 1]),
)
h = (
np.sin((phi2 - phi1) / 2) ** 2
+ np.cos(phi1) * np.cos(phi2) * np.sin(dlng / 2) ** 2
)
table["distance_km"] = 2 * 6371 * np.arcsin(np.sqrt(h))
return fastmob.Trips(table)
half = set(
np.random.default_rng(0).permutation(us["uid"].unique())[: us["uid"].nunique() // 2]
)
observed_flows = fastmob.FlowDataFrame(observed, locations=catalogue)
predicted_flows = fastmob.FlowDataFrame(predicted, locations=catalogue)
cpc = evaluation.common_part_of_commuters(observed_flows, predicted_flows)
cpl = evaluation.common_part_of_links(observed_flows, predicted_flows)
trips_a = with_distance(
fastmob.trajectory_to_trips(us[us["uid"].isin(half)], resolution=7)
)
trips_b = with_distance(
fastmob.trajectory_to_trips(us[~us["uid"].isin(half)], resolution=7)
)
cpcd = evaluation.common_part_of_commuters_distance(trips_a, trips_b)
cpcd_same = evaluation.common_part_of_commuters_distance(trips_a, trips_a)
cpcd_method = trips_a.common_part_of_commuters_distance(trips_b)
print(f"observed vs gravity: CPC {cpc:.2f}, CPL {cpl:.2f}")
print(
f"CPCD between two halves of the users: {cpcd:.2f} (as a Trips method: {cpcd_method:.2f}); of a set of trips with itself: {cpcd_same:.2f}"
)
observed vs gravity: CPC 0.68, CPL 0.34
CPCD between two halves of the users: 0.86 (as a Trips method: 0.86); of a set of trips with itself: 1.00
About 68% of the trips follow links present in both tables (CPC 0.68), but only 34% of the links are shared (CPL 0.34): the model puts flow on many links where little or none was observed. The CPCD of two halves of the users is 0.86 (similar trip-distance distributions; the Trips method gives the same value), and it is exactly 1.00 when a set of trips is compared with itself, as it should be.
Compare by group
compare_to applies a metric to a column of two dataframes, overall and optionally per group, and returns a ComparisonResult with overall (the mean over groups) and the per-group values. The default metric is the Wasserstein distance; pass metric= for another, such as jensen_shannon_divergence. Here, flow values grouped by origin zone.
result = evaluation.compare_to(
pairs.rename(columns={"flow_obs": "flow"})[["origin", "destination", "flow"]],
pairs.rename(columns={"flow_pred": "flow"})[["origin", "destination", "flow"]],
"flow",
group_col="origin",
)
print(type(result).__name__, "- overall", round(result.overall, 2))
worst = sorted(result.groups, key=lambda g: -g[1])[:3]
print(
"zones with the largest deviation:",
[(zone, round(value, 1)) for zone, value in worst],
)
ComparisonResult - overall 1.04
zones with the largest deviation: [('53', 13.7), ('40', 9.1), ('63', 8.8)]
The overall value, 1.04, is the mean over origin zones of the Wasserstein distance between observed and predicted flow values; the zones with the largest deviation are 53 (13.7), 40 (9.1) and 63 (8.8). It points to the origins where the model fits worst.
Compare when things happen: spatio-temporal distance
stvd_emd compares two spatio-temporal visitation distributions (how many people are at each place in each 10-minute bin of the day): an entropy-regularised earth mover's distance between them, so moving mass to a nearby place or time costs less than to a distant one. The cost is in metres, and one minute of time counts as alpha = 10 metres. We compare the occupancy of two random halves of the users (each check-in assumed to last 30 minutes), on an H3 grid of resolution 6 (about 36 km²). It is the slowest metric here.
stays = fastmob.Staypoints(
pd.DataFrame(
{
"uid": us["uid"],
"lat": us["lat"],
"lng": us["lng"],
"started_at": us["datetime"],
"finished_at": us["datetime"] + pd.Timedelta(minutes=30),
}
)
)
places, stays_at = stays.generate_global_locations(h3_resolution=6)
table = stays_at.to_pandas()
occupancy_a = fastmob.Staypoints(table[table["uid"].isin(half)]).build_stvd(places)
occupancy_b = fastmob.Staypoints(table[~table["uid"].isin(half)]).build_stvd(places)
emd = evaluation.stvd_emd(occupancy_a, occupancy_b)
print(
f"{len(occupancy_a)} and {len(occupancy_b)} place-time cells; spatio-temporal distance {emd:,.0f}"
)
The two halves of the users have 72,922 and 72,967 place-time cells and a spatio-temporal distance of about 23,400 (metres of transport cost, with one minute counting as 10 m). It is large because each half has only about 150 users spread across the continent; two halves of a much larger sample would be closer. Use it to rank alternatives (a model against data, one preprocessing against another), not as an absolute measure.
Interpretation and limits
- No single metric is "the" accuracy. Error metrics (
rmse,r_squared) are dominated by the largest flows, shape metrics (jensen_shannon_divergence,wasserstein_distance) by the many small ones, and CPC by how much traffic follows shared links. Report several, chosen for your question. - The gravity model is fitted and judged on the same flows, so these numbers say how well it describes the data, not how well it predicts new data. Fit on one part of the data and test on another for that.
- Zero flows matter: a zone pair with predicted flow but no observed flow is an error, and the shape metrics need a small constant to avoid taking the logarithm of zero.
- CPC and CPL depend on the zone size: the same data compared at finer zones will agree less.
- The spatio-temporal distance depends on the grid, the stay duration and its parameters; its absolute value is not meaningful on its own, only differences between comparisons are.
Next steps
- Generate trajectories to compare against: Generate mobility with models.
- Fit laws to your data: Compare your data to mobility laws.
Related API
mse ·
rmse ·
nrmse ·
r_squared ·
max_error ·
kullback_leibler_divergence ·
information_gain ·
jensen_shannon_divergence ·
wasserstein_distance ·
common_part_of_commuters ·
common_part_of_links ·
common_part_of_commuters_distance ·
compare_to ·
stvd_emd