Experimental API reference

These modules diagnose dependence and validate collection designs. Their interfaces may change before the stable inference API does.

Spatial and temporal dependence diagnostics for geoestimate.

These experimental functions measure within-cluster correlation and how fast it decays along geographic distance or elapsed time. They describe dependence; they do not replace the standard-error estimators in geoestimate.inference. other pairwise-distance matrix. They are diagnostics: they describe the dependence structure (and how much information itineraries cost), they do not replace the standard-error estimators in inference.py.

The core is axis-agnostic. empirical_variogram, morans_i, and effective_n all take a precomputed n×n distance matrix, so the same code serves the spatial axis (haversine_matrix), the temporal axis (time_gap_matrix), and any future axis (e.g. same/different enumerator).

References

Griffith, D.A. (2005). Effective geographic sample size in the presence

of spatial autocorrelation.

Watson, P.A. (2021). A note on the variogram-based effective sample size.

10. Applied Statistics.

class geoestimate.spatial.DependenceDiagnostics(within_between_ratio=nan, morans_i_space=nan, morans_i_space_p=nan, variogram_range_m=nan, spatial_corr_ratio=nan, n_eff_space=nan, morans_i_time=nan, morans_i_time_p=nan, variogram_range_s=nan, temporal_corr_ratio=nan, n_eff_time=nan)

Experimental spatial and temporal dependence diagnostics.

Parameters:

• within_between_ratio (float)

• morans_i_space (float)

• morans_i_space_p (float)

• variogram_range_m (float)

• spatial_corr_ratio (float)

• n_eff_space (float)

• morans_i_time (float)

• morans_i_time_p (float)

• variogram_range_s (float)

• temporal_corr_ratio (float)

• n_eff_time (float)

geoestimate.spatial.haversine_matrix(lon, lat)

Pairwise great-circle distances in meters between lon/lat points.

Parameters:

• lon (ndarray) – Longitudes in decimal degrees.

• lat (ndarray) – Latitudes in decimal degrees.

Returns:

n×n symmetric matrix of distances in meters (zero diagonal).

Raises:

ValueError – If coordinates are misaligned, nonfinite, or outside valid longitude and latitude bounds.

Return type:

ndarray

geoestimate.spatial.time_gap_matrix(timestamps)

Pairwise absolute time differences in seconds.

Parameters:

timestamps (ndarray) – datetime64 / Timestamps / ISO strings, or numeric epoch seconds.

Returns:

n×n symmetric matrix of absolute time gaps in seconds (zero diagonal).

Return type:

ndarray

geoestimate.spatial.empirical_variogram(values, dist, n_bins=15, max_dist=None)

Classical (Matheron) empirical semivariogram over distance bins.

Semivariance per bin is the mean of 0.5 * (z_i - z_j)**2 over unique pairs whose distance falls in the bin.

Parameters:

• values (ndarray) – Length-n outcome array.

• dist (ndarray) – n×n pairwise-distance matrix (any axis).

• n_bins (int) – Number of distance bins.

• max_dist (float | None) – Maximum distance to include. Defaults to half the maximum pairwise distance (the conventional cutoff).

Returns:

(lag_centers, semivariance, pair_counts) for non-empty bins only.

Raises:

ValueError – If inputs are malformed or outside their domains.

Return type:

tuple[ndarray, ndarray, ndarray]

geoestimate.spatial.fit_variogram(lags, gamma, counts, model='exponential')

Weighted least-squares fit of an exponential variogram.

The fitted correlation function is rho(h) = ((C1-C0)/C1) * exp(-h/r), so r is the e-folding scale (correlation falls to 1/e at h = r); the conventional “effective range” is ~3r.

Parameters:

• lags (ndarray) – Bin-center distances from empirical_variogram.

• gamma (ndarray) – Semivariances from empirical_variogram.

• counts (ndarray) – Pair counts per bin (used as fit weights).

• model (str) – Only "exponential" is supported.

Returns:

(nugget C0, sill C1, range r). NaNs if the fit cannot be performed.

Raises:

ValueError – If model names anything but "exponential".

Return type:

tuple[float, float, float]

geoestimate.spatial.morans_i(values, dist, cutoff, n_perm=999, seed=0)

Global Moran’s I with binary distance-cutoff weights.

Weight w_ij = 1 if 0 < dist_ij <= cutoff else 0. The p-value is a two-sided permutation test on the observed statistic.

Parameters:

• values (ndarray) – Length-n outcome array.

• dist (ndarray) – n×n pairwise-distance matrix.

• cutoff (float) – Neighbor distance threshold (same units as dist).

• n_perm (int) – Number of permutations for the p-value.

• seed (int) – RNG seed for the permutation test.

Returns:

(I, p_value). NaNs if undefined (e.g. no neighbor pairs).

Raises:

• TypeError – If seed is not an integer.

• ValueError – If inputs are malformed or outside their domains.

Return type:

tuple[float, float]

geoestimate.spatial.effective_n(values, dist, nugget, sill, range_)

Variogram-based effective sample size under autocorrelation.

n_eff = n / (1 + (1/n) * sum_{i!=j} rho(d_ij)) with the exponential correlation rho(d) = ((sill - nugget)/sill) * exp(-d/range_) (Griffith 2005; Watson 2021). Reduces to n when there is no autocorrelation.

Parameters:

• values (ndarray) – The observations, used only for their count.

• dist (ndarray) – Square pairwise distance matrix over those observations.

• nugget (float) – Variogram nugget C0.

• sill (float) – Variogram sill C1.

• range – Variogram e-folding range r.

• range_ (float)

Returns:

Effective sample size (<= n), or n if the fit was degenerate.

Raises:

ValueError – If the distance matrix or variogram parameters are outside their supported domains.

Return type:

float

geoestimate.spatial.within_between_contrast(values, labels)

Mean semivariance of same-cluster vs different-cluster pairs.

This is the “does the prior point predict the next one, versus a point elsewhere?” contrast. ratio = within / between near 1 means itineraries look like representative subsamples (dispersed); ratio well below 1 means within-itinerary pairs are much more alike (compact routes), which inflates the design effect.

Parameters:

• values (ndarray) – One value per observation.

• labels (ndarray) – Cluster label per observation, aligned to values.

Returns:

{"within", "between", "ratio"} (semivariances; NaN where undefined).

Raises:

ValueError – If values and labels are malformed or misaligned.

Return type:

dict[str, float]

geoestimate.spatial.dependence_diagnostics(values, labels, *, lon=None, lat=None, timestamps=None, max_points=2500, seed=0)

Compute experimental within-cluster dependence diagnostics.

Parameters:

• values (ndarray) – One finite outcome per observation.

• labels (ndarray) – Cluster labels aligned with values.

• lon (ndarray | None) – Optional longitudes in decimal degrees.

• lat (ndarray | None) – Optional latitudes in decimal degrees.

• timestamps (ndarray | None) – Optional datetimes or epoch seconds.

• max_points (int) – Maximum observations used for pairwise diagnostics.

• seed (int) – Random seed used for subsampling and permutations.

Returns:

Spatial and temporal descriptive diagnostics.

Raises:

• TypeError – If seed is not an integer.

• ValueError – If arrays are misaligned, incomplete, or outside their documented domains.

Return type:

DependenceDiagnostics

Space-time Monte Carlo harness for stress-testing itinerary-based collection.

The point of this module is limiting-case reasoning. We simulate a known space-time data-generating process (DGP) with no measurement noise, run real collection pipelines (sampling + routing under a fixed per-route time budget) on top, push the assumptions to their extremes, and check what happens to the bias, standard-error calibration, and confidence-interval coverage that geoestimate.Sample produces.

The analytically-predictable limits (the “truth table”):

Spatial correlation length range_s_m vs. the route’s inter-point spacing

• 0 (white noise): information ∝ point count → COMPACT routing is BLUE; ICC≈0, deff≈1, naive SE correct.

• ∞ (field constant in space): n_eff_space → 1; extra spatial points carry no information; only the time axis can buy precision.

Temporal structure (diurnal_amp, range_t_min, start-time policy)

• none: start-time policy is irrelevant, no temporal bias.

• strong diurnal + synchronized starts: β̂ is biased toward the sampled time-of-day by a fixed amount; more points do NOT fix it; staggering start times across the day does.

Key separations the extremes make obvious:

• Correlation moves VARIANCE / SE / coverage, never the bias of the mean.

• Bias comes from SELECTION: spatial coverage gaps and non-representative time windows.

• The coverage cliff (naive SE under compact routing) is an INTERIOR phenomenon: benign at both correlation extremes, worst in the middle.

Everything works in lon/lat with great-circle distances so it exercises the same code path as real data.

class geoestimate.simulate.SimConfig(grid_n=18, extent_deg=0.05, range_s_m=800.0, sd_s=1.0, base_logit=-0.85, day_min=600.0, range_t_min=60.0, diurnal_amp=0.0, sd_t=0.0, time_grid_n=40, n_itineraries=8, shift_min=120.0, speed_m_per_min=80.0, dwell_min=2.0, n_sims=200, seed=12345)

Parameters of the space-time DGP and the field operation.

Defaults describe a ~5.5 km square region surveyed in shifts of a few hours. Push any parameter to an extreme to read off a limiting case.

Parameters:

• grid_n (int)

• extent_deg (float)

• range_s_m (float)

• sd_s (float)

• base_logit (float)

• day_min (float)

• range_t_min (float)

• diurnal_amp (float)

• sd_t (float)

• time_grid_n (int)

• n_itineraries (int)

• shift_min (float)

• speed_m_per_min (float)

• dwell_min (float)

• n_sims (int)

• seed (int)

grid_side()

Number of candidate points per grid side (grid is this squared).

Return type:

int

class geoestimate.simulate.Population(lon, lat, dist_m, g_s, t_grid, temporal, base_logit, beta_true)

One realized space-time population (no measurement noise).

Parameters:

• lon (ndarray)

• lat (ndarray)

• dist_m (ndarray)

• g_s (ndarray)

• t_grid (ndarray)

• temporal (ndarray)

• base_logit (float)

• beta_true (float)

p_at(idx, times_min)

True ratio p(s_i, t_i) at points idx visited at times_min.

Parameters:

• idx (ndarray)

• times_min (ndarray)

Return type:

ndarray

class geoestimate.simulate.PopulationFactory(cfg, lon=None, lat=None)

Precomputes the geometry + covariance factors, then draws populations.

Factoring once per configuration keeps the Monte Carlo loop cheap: each draw is just two matrix-vector products.

Parameters:

• cfg (SimConfig)

• lon (ndarray | None)

• lat (ndarray | None)

draw(rng)

Draw one space-time population and its true space-time-average ratio.

Parameters:

rng (Generator)

Return type:

Population

class geoestimate.simulate.Pipeline(name, routing='compact', staggered_starts=True)

A named collection strategy.

Parameters:

• routing (str) – How the next point is chosen. “compact” selects the nearest unvisited, “dispersed” (farthest from already-visited, i.e. maximin), “systematic” (a fixed space-filling order from a random offset), or “srs” (uniformly random next).

• staggered_starts (bool) – If True, route start times are spread across the day; if False, all routes start at the same time-of-day (synchronized).

• name (str)

geoestimate.simulate.collect(pop, pipe, cfg, rng)

Run all K routes of a pipeline and return an annotated-frame DataFrame.

Parameters:

• pop (Population)

• pipe (Pipeline)

• cfg (SimConfig)

• rng (Generator)

Return type:

DataFrame

class geoestimate.simulate.PipelineResult(name, se_method, n_sims, mean_n, mean_dist_km, bias, true_sd, mean_se, se_sd_ratio, coverage, mean_n_eff_space, mean_within_between)

Monte Carlo metrics for one (pipeline, SE method) combination.

bias and true_sd describe the sampling distribution of beta-hat across sims; se_sd_ratio (mean SE / true SD) and coverage describe whether the reported uncertainty is honest.

Parameters:

• name (str)

• se_method (str)

• n_sims (int)

• mean_n (float)

• mean_dist_km (float)

• bias (float)

• true_sd (float)

• mean_se (float)

• se_sd_ratio (float)

• coverage (float)

• mean_n_eff_space (float)

• mean_within_between (float)

as_row()

Return a flat, rounded dict of metrics for tabular display.

Return type:

dict[str, float | str | int]

geoestimate.simulate.run_pipeline(factory, pipe, cfg, se_method='auto', ci_level=0.95, spatial_diag=False)

Monte Carlo a single synthetic pipeline; return aggregated metrics.

Parameters:

• factory (PopulationFactory) – Precomputed population factory (geometry + covariance).

• pipe (Pipeline) – The collection strategy to simulate (synthetic routing).

• cfg (SimConfig) – Simulation configuration (DGP + field operation + Monte Carlo).

• se_method (str) – Which standard error drives the CI: “auto”, “naive”, “cluster” (analytic robust with t_{G-1}), or “boot” (pairs bootstrap).

• ci_level (float) – Nominal confidence level.

• spatial_diag (bool) – If True, also collect the spatial dependence diagnostics.

Returns:

A PipelineResult over cfg.n_sims simulations.

Raises:

ValueError – If the method or confidence level is unsupported.

Return type:

PipelineResult

geoestimate.simulate.evaluate_scene(factory, sample_idx, itinerary_id, time_of_day_min, timestamp_s, cfg, se_method='auto', ci_level=0.95, spatial_diag=True, label='scene')

Design-conditional Monte Carlo of a realized survey against the city mean.

factory carries the field on the whole city universe; the survey is a fixed subsample of it (sample_idx) with a fixed itinerary partition and visit times (you ran geo_sampling + the allocator once). Each simulation redraws the outcome field, estimates the ratio from the observed sample, and checks whether the CI covers the universe space-time mean. Because the sample is a strict subset of the universe, β̂ carries genuine spatial sampling error. The cluster SE must capture that error, so this is a real test of whether estimate’s SE/CI is honest for that design.

Parameters:

• factory (PopulationFactory) – PopulationFactory(cfg, lon, lat) built on ALL city points.

• sample_idx (ndarray) – Indices (into the factory universe) of the observed frames.

• itinerary_id (ndarray) – Per-frame cluster labels, aligned to sample_idx.

• time_of_day_min (ndarray) – Per-frame time-of-day in [0, day_min] minutes. This drives the diurnal/temporal field component.

• timestamp_s (ndarray) – Per-frame absolute timestamp in seconds across the whole operation, used by the optional temporal diagnostics.

• cfg (SimConfig) – Simulation configuration (DGP + Monte Carlo).

• se_method (str) – SE estimator to use, or “auto” to take the design’s.

• ci_level (float) – Coverage the interval claims.

• spatial_diag (bool) – Whether to accumulate the spatial dependence diagnostics.

• label (str) – Display name for the result row.

Returns:

A PipelineResult over cfg.n_sims field redraws.

Raises:

ValueError – If the method, confidence level, or scene arrays are invalid.

Return type:

PipelineResult

geoestimate.simulate.run_experiment(cfg, pipelines, se_method='auto', ci_level=0.95, spatial_diag=False)

Run every pipeline against one shared population factory (one DGP).

Parameters:

• cfg (SimConfig)

• pipelines (list[Pipeline])

• se_method (str)

• ci_level (float)

• spatial_diag (bool)

Return type:

list[PipelineResult]

geoestimate.simulate.results_table(results)

Render results as a fixed-width table.

Parameters:

results (list[PipelineResult])

Return type:

str

geoestimate.simulate.default_pipelines()

The four reference routing strategies, all with staggered start times.

Return type:

list[Pipeline]

geoestimate.simulate.main()

Reproduce the headline experiments and print them as tables.

Run with python -m geoestimate.simulate.

Return type:

None

Bridge from the real geosensing pipeline to a geoestimate simulation scene.

This turns the actual geo_sampling to allocator output into a fixed scene with point coordinates, itinerary membership, and visit times. The scene geoestimate.simulate.evaluate_scene can validate a DGP against, and that mirrors what the annotated frames look like in production.

allocator and geo_sampling are optional; install them with pip install geoestimate[pipeline] (or uv pip install -e ../allocator ../geo_sampling for local checkouts). They are imported lazily so core geoestimate keeps no heavy dependencies.

class geoestimate.pipeline.Scene(longitude, latitude, itinerary_id, timestamp_s, time_of_day_min)

A realized field operation: where, in which itinerary, and when.

All arrays are per annotated frame and aligned by position.

Parameters:

• longitude (ndarray)

• latitude (ndarray)

• itinerary_id (ndarray)

• timestamp_s (ndarray)

• time_of_day_min (ndarray)

to_frame()

Annotated-frame layout (minus the outcome columns).

Return type:

DataFrame

property n_itineraries: int

Count the distinct itineraries the frames belong to.

Returns:

The number of unique itinerary ids.

property day_span: float

Number of days the operation spans (from absolute timestamps).

geoestimate.pipeline.points_from_roads(roads, per_segment=1)

Point locations along road segments, from the road-segment schema.

Accepts a DataFrame or a CSV path with the geo_sampling / allocator columns start_lat, start_long, end_lat, end_long (and returns longitude/latitude). If the frame already has longitude / latitude it is returned unchanged.

per_segment interpolates that many evenly-spaced points along each segment (1 = midpoint). Densifying gives a larger candidate universe and hence a small sampling fraction for realistic validation.

Parameters:

• roads (DataFrame | str)

• per_segment (int)

Return type:

DataFrame

geoestimate.pipeline.build_itineraries(points, method='random_partition', n_itineraries=None, max_distance=None, seed=None)

Partition points into itineraries with allocator (offline haversine).

Parameters:

• points (DataFrame) – Table containing longitude and latitude columns.

• method (str) – Allocator partitioning method.

• n_itineraries (int | None) – Requested itinerary count, when supported by the method.

• max_distance (float | None) – Optional maximum route distance in meters.

• seed (int | None) – Optional allocator random seed.

Returns:

The allocator data table and routes as point-index lists in visit order.

Raises:

• ImportError – If the optional allocator package is unavailable.

• TypeError – If seed is not an integer or None.

• ValueError – If points or options are outside their supported domains.

• RuntimeError – If allocator does not return a complete partition.

Return type:

tuple[DataFrame, list[list[int]]]

geoestimate.pipeline.assign_visit_times(points, routes, *, days=14, shifts_per_day=1, speed_m_per_min=80.0, dwell_min=2.0, day_minutes=600.0, stagger_starts=True, seed=None)

Spread itineraries over a multi-day operation and time every frame.

Each itinerary is a single shift assigned to one (day, slot). Within a shift, visit times accumulate along the route as travel (haversine metres / speed_m_per_min) plus a per-point dwell_min. Shifts are spread round-robin across days × shifts_per_day so timestamps span the whole operation while each shift’s frames stay close in time-of-day.

Returns a Scene whose time_of_day_min drives the diurnal field and whose timestamp_s is the absolute time fed to estimate’s temporal diagnostic.

Parameters:

• points (DataFrame)

• routes (list[list[int]])

• days (int)

• shifts_per_day (int)

• speed_m_per_min (float)

• dwell_min (float)

• day_minutes (float)

• stagger_starts (bool)

• seed (int | None)

Return type:

Scene

geoestimate.pipeline.make_scene(roads, method='random_partition', n_itineraries=200, seed=0, **time_kwargs)

Convenience: roads CSV/frame → points → itineraries → timed Scene.

Parameters:

• roads (DataFrame | str)

• method (str)

• n_itineraries (int)

• seed (int | None)

• time_kwargs (object)

Return type:

Scene

geoestimate.pipeline.subsample_scene(universe, n_sample, method='kmeans_tsp', n_itineraries=80, seed=0, stagger_starts=True, days=14, shifts_per_day=None)

Sample a survey out of a city universe and route it into itineraries.

Treats universe (all candidate road segments for a city) as the population, draws n_sample of them (SRS as a stand-in for live geo_sampling), and routes the sample with the allocator. Returns (sample_idx, scene) where sample_idx indexes the universe (so the field can be drawn on the whole city and the sample scored against the city mean in geoestimate.simulate.evaluate_scene).

Parameters:

• universe (DataFrame | str)

• n_sample (int)

• method (str)

• n_itineraries (int)

• seed (int)

• stagger_starts (bool)

• days (int)

• shifts_per_day (int | None)

Return type:

tuple[ndarray, Scene]