How the fits work

Question. How does Fastkit-Mobility estimate a distance law or a gravity model, how sure can I be of the result, and what must I report so that a reader can check it?

This page explains the estimators, their limits and a reporting checklist. To run them, open Compare to mobility laws for distributions, Validate against census flows for flows, or Validate mobility with GeoLife for the whole cycle on GPS tracks: trip-length law, gravity distance decay and a gravity against radiation hold-out.

Distributions of distances

fit_distance_distribution fits one family to positive values (jump lengths, trip lengths, radii of gyration) by maximum likelihood. compare_distance_distributions fits several and ranks them.

Step	What is done, and why it matters
Range	Values outside `[xmin, xmax]`, and zeros, are left out and counted in `n_excluded`. Default: smallest positive value to largest. Estimates depend on the range. Zeros are common (repeated check-ins at one place) and have no place in a power law. Raise `xmin` to ignore GPS noise, or let `select_xmin` pick it by the KS distance (Clauset et al., 2009). Set `xmax` to the largest distance that is possible in the study area (the diagonal of its bounding box, or half the Earth's circumference for global data): if it is left at the sample maximum, the sample maximum is the end of the distribution and a cutoff can hide behind it, so the library warns.
Grouping	By default the values are grouped into 200 geometric bins on the range, and the likelihood is computed from the bin counts, so a million values cost no more than a thousand. `likelihood="exact"` uses the individual values instead (no binning approximation, but slow beyond about 10^5 values).
Likelihood	Each family is normalised on the same range, and the multinomial log-likelihood of the bin counts is maximised. Families are comparable: they explain the same data on the same range.
Uncertainty	The bin counts are resampled (bootstrap, 100 replicates by default) and the model refitted; the report has the standard error and the percentile interval of each parameter. Values from the same person are not independent: pass `groups=` (for example the user of each jump) to resample whole people, which gives wider and more honest intervals.
Fit quality	The KS distance between the empirical and the fitted cumulative distribution; with `n_gof` above 0, a parametric-bootstrap p-value. See the caveat below.
Comparison	AIC and BIC, their difference to the best model, and the Akaike weight; the likelihood-ratio test for nested families and the Vuong test for the others.
Search range	Each parameter is searched within a range (`bounds=` changes it). A fit that ends on the edge of its range is flagged.

The five families are the power law x^-beta, the shifted power law (x + r0)^-beta, the exponential exp(-x / kappa), the log-normal, and the truncated power law (x + r0)^-beta exp(-x / kappa) of González et al. (2008). The shifted and the truncated power law differ only by the cutoff kappa, so comparing them shows whether the data support a cutoff at all.

Reading the result

A parameter at its bound. Each parameter is searched within a range. If an estimate lands on the edge, at_bound lists it and the data do not identify it. Some edges are legitimate limits of a family, not a distortion (DistributionFit.limits): r0 at its lower edge turns the shifted power law into a power law and the truncated power law into a power law with an exponential cutoff, and kappa at its upper edge (1,000 times xmax by default) turns the truncated power law into the shifted one. Every other edge, for example a shifted power law with beta at its upper limit, is degenerate (DistributionFit.degenerate). For the truncated power law, a kappa at its upper bound means that, below xmax, there is no cutoff in the data: the model is a shifted power law, and a fitted kappa of 10^8 km is not evidence for a cutoff. Its standard error and interval are reported as NaN: say that the cutoff is not identified. fit.profile("kappa") then gives what the data do say, a profile-likelihood interval (each end found by bisection, so it does not depend on a grid) whose open end reads "kappa is at least X km"; its nuisance_at_bound lists other parameters that sit on an edge along the profile, which means that a degenerate shape helps to decide the interval; use bounds= to show that the conclusion does not depend on the search range.
Independence. The intervals and the p-value assume that the values are independent. Jumps from the same person are not, so intervals computed without groups= understate the uncertainty about people in general (they describe the fitting error only).
Goodness of fit. With thousands of observations every simple model is rejected, because the test detects tiny departures. Report the KS distance and the AIC difference, and look at the fit on the data, instead of reading the p-value alone. A small p-value does not mean the model is useless, and a large one does not prove it right.
Does the data need a cutoff? likelihood_ratio_test(shifted, truncated) asks it directly. The cutoff is a limit of the model (kappa to infinity), so the p-value uses the boundary mixture of chi-squared distributions, not the plain one, and that asymptotic p-value is too small when xmax clips the data or the test is far in the tail; n_bootstrap= gives a parametric-bootstrap p-value, which is the one to report. The same function tests powerlaw against shifted_powerlaw.
Families that are not nested. vuong_test(a, b) asks whether two families (for example a log-normal and a truncated power law) fit the data differently; a small statistic means the data cannot choose, whatever the AIC difference. Like the other tests it treats the values as independent.
Several families. An AIC difference below about 2 means the data cannot choose. The Akaike weight is a relative support among the families you list, not a probability that the family is true. The likelihood is on binned counts: use differences between families, not the absolute value. Read the at_bound column first: a family with a parameter on the edge of its range (for example a shifted power law whose beta is at its upper limit) is a degenerate fit, and its AIC rank means little. compare_distance_distributions therefore lists degenerate families (not the legitimate limits above) last (degenerate column), measures delta_aic from the best family that is not degenerate and leaves them out of the Akaike weights; pass degenerate_last=False for the plain ranking. The AIC gaps and weights also treat the values as independent, so they overstate the evidence when many values come from the same people.
Which estimate to report. The histogram grid search and the maximum-likelihood fit give different exponents on the same data (for the check-in jumps 1.61 against 1.76, because their ranges and estimators differ). Report the maximum-likelihood one, with its range and interval.
The range. A single fit is conditional on xmin. select_xmin fits a grid of candidates and picks the one with the smallest KS distance; it returns the whole table, so report how the estimate moves with xmin. The intervals of the final fit do not include the uncertainty about the selected xmin.
The earlier fit. fit_values_to_truncated_powerlaw fits a log histogram by a grid search and returns no uncertainty. It is kept for the plotting helpers; use fit_distance_distribution for a result that you report.

Gravity models

Gravity.fit fits flow ∝ origin^a * destination^b * f(distance) by Poisson regression with a log link (Flowerdew and Murray, 1982), the standard estimator for flow counts (Santos Silva and Tenreyro, 2006, discuss its use for gravity models).

Output	Meaning
`fit_summary()`	Estimate, standard error, 95% confidence interval and p-value of each parameter. The distance row is `distance_coefficient`, the exponent in `flow ∝ distance^gamma` (negative), or `distance_decay` for the exponential form, in km⁻¹.
`fit_statistics()`	Number of flows, log-likelihood, deviance, pseudo R-squared (one minus the deviance over the null deviance) and the dispersion.

Robust errors. Flow data are almost always overdispersed: their variance is larger than their mean, so the standard errors of a Poisson model are too small. The dispersion in fit_statistics() (Pearson chi-squared per degree of freedom) is 1 for a good Poisson fit and well above 1 for real flows (about 16 for the check-in flows in the notebook). On a table of expanded sample counts it also reflects the size of the weights, so do not read it as behavioural overdispersion there. Gravity.fit therefore reports heteroskedasticity-robust (sandwich) errors by default; the estimates do not change. Pass robust=False for the model-based errors, or choose the estimator with cov_type ("HC0", "HC1" or "nonrobust"). HC1 multiplies HC0 by n / (n - k) (the small-sample factor of Stata's robust option and R's sandwich::vcovHC(type = "HC1")), with the origin levels counted in k for the singly constrained model; use it to match an analysis made with those tools. fit_statistics() reports the estimator used.
Singly or globally constrained. The singly constrained model has one free level for each origin, so it matches the outflow of every origin and the distance exponent is estimated from how flows spread among destinations. It is the right choice when the origin totals are observed (a census, or the origins of a sample). The globally constrained model has a single scale and fits only if the total flow of each origin follows the model.
Sizes of origin and destination. By default one catalogue column (relevance_column) is the size of both sides. Pass origin_relevance_column (for example residents or departures) and destination_relevance_column (for example jobs or arrivals) to give them separate sizes. The origin size has an exponent only in the globally constrained model; in the singly constrained one the origin is absorbed in its own level, so Gravity.fit raises an error if you name an origin column there.
Zeros. A table of observed flows (such as the output of trajectory_to_od) lists only the pairs that have flow. Fitting those alone leaves out the long pairs that nobody travelled, and flattens the distance decay (for the check-in flows in the notebook, -1.21 without them and -1.40 with them). Gravity.fit therefore adds every missing pair of the catalogue as a flow of 0 (include_zeros=True), for origins that have at least one flow. The singly constrained fit keeps one row per pair, so 1,000 locations take about 400 MB and 10 seconds; the origin levels are solved without building dummy columns. This is the point of the Poisson estimator that Santos Silva and Tenreyro (2006) recommend.
Unweighted samples. A bias that depends only on the origin (more people sampled in some home zones) does not bias the distance exponent of a singly constrained fit when the model holds, but the weights do change the estimate and its precision in a finite sample, and they change the flow levels. Weight the sample before comparing absolute flows (see the notebook).
Pseudo R-squared. It compares the deviance with that of a model with only a constant. In the singly constrained model most of the explained deviance is the size of the origins, so a high value does not show that the distance law holds: look at the distance coefficient and its interval.
Relevance. If the relevance of a destination is made of the same flows (for example check-ins counted at the destination), the destination exponent is partly circular. Prefer an independent size such as population or jobs.
Zones. The exponent depends on the zone size and on the range of distances in the data. Report both.

Comparing flows

compare_flow_models scores model flows against observed flows over a common set of origin-destination pairs.

Score	Meaning
`cpc`	Common part of commuters: `2 * sum(min(observed, model)) / (sum(observed) + sum(model))`. 1 is identical.
`cpl`	Common part of links: the same for the sets of pairs that have a flow, whatever its size.
`r_squared`, `spearman`	Squared Pearson and rank correlation of the pair flows.
`poisson_deviance`, `pseudo_r_squared`	Fit of the model flows as the mean of Poisson counts, against a model that predicts the average flow everywhere.

All models are scored on the same set of pairs, by default every pair of the observed table's locations, so the rows can be compared and a model's score does not depend on which other models are scored with it (report n_pairs). With pairs="union" the set is the pairs that have flow in any table, and then adding a model changes the others' scores. All scores depend on the zone size. CPC is dominated by the largest flows, so a high CPC can hide a poor fit to the many small ones: look at the rank correlation too. CPL is only informative for integer (sampled) flows: expected flows are positive on every pair, which fixes their CPL. A model that was fitted to the observed flows is scored in-sample, so compare it with a parameter-free model (radiation), or score it on flows that were held out, as the notebook does.

What to report

The data, the filtering, and the number of values or flows used.
For distributions: the family, the range [xmin, xmax] and n_excluded, the estimates with their intervals, the number of bootstrap replicates and whether people were resampled (groups=), the parameters at a bound, the KS distance (and the p-value with its replicates), the search bounds if you changed them, and the comparison with the alternatives (AIC and weights, the likelihood-ratio or Vuong test, and which fits were degenerate).
For gravity fits: the deterrence function, the constraint type, whether zero pairs were included, the estimates with intervals, whether the errors are robust, the dispersion, the zone size and the source of the relevance.
For validation against a reference table: the weighting, the scores and the zone size, with the same scores for a trivial baseline (for example the unweighted sample).
The version of Fastkit-Mobility and the random seeds.

References

González, M. C., Hidalgo, C. A. & Barabási, A.-L. (2008). Understanding individual human mobility patterns. Nature 453, 779-782.
Clauset, A., Shalizi, C. R. & Newman, M. E. J. (2009). Power-law distributions in empirical data. SIAM Review 51, 661-703.
Flowerdew, R. & Murray, A. (1982). A method of fitting the gravity model based on the Poisson distribution. Journal of Regional Science 22, 191-202.
Santos Silva, J. M. C. & Tenreyro, S. (2006). The log of gravity. The Review of Economics and Statistics 88, 641-658.
Chernoff, H. (1954). On the distribution of the likelihood ratio. Annals of Mathematical Statistics 25, 573-578.
Vuong, Q. H. (1989). Likelihood ratio tests for model selection and non-nested hypotheses. Econometrica 57, 307-333.
Simini, F., González, M. C., Maritan, A. & Barabási, A.-L. (2012). A universal model for mobility and migration patterns. Nature 484, 96-100.

Next steps

Fit and compare the distributions of your own jump lengths: Compare to mobility laws.
Check flows against a reference table: Validate against census flows.
The whole cycle on GPS tracks: Validate mobility with GeoLife.

fit_distance_distribution · compare_distance_distributions · DistributionFit · select_xmin · likelihood_ratio_test · vuong_test · Gravity.fit · compare_flow_models