bayesian-pv-census

Turn a detector's raw count into a census with credible intervals, then audit a register against it.

A detector that finds objects in imagery misses some and invents others, at rates that vary from place to place. Its raw total is therefore not a measurement, and comparing it directly to an official register says as much about the detector as about the register. This package closes that gap: given a validation sample, it turns the raw total into a posterior over the true quantity, and turns that posterior into a verdict on whatever the register reports.

Nothing here knows about photovoltaics, France, or geometry. A unit is anything with a raw total and a validation sample — a department, a grid cell, a utility service area. raw is whatever you chose to count: installed capacity, number of installations, roof area.

This is the statistical core of Nationally Consistent, Locally Incomplete: A Bayesian Remote-Sensing Audit of Rooftop Photovoltaic Registries (Kasmi et al., 2026, pending peer review) , extracted as a library. Every public function corresponds to a named component of the paper.

Install

pip install bayesian-pv-census            # numpy, scipy, pandas
pip install 'bayesian-pv-census[figures]' # adds matplotlib

One minute

from bayesian_pv_census import UnitRecord, correct_unit, evaluate

unit = UnitRecord(
    unit_id="dept_86",
    raw=26_800,                          # what the detector found, any consistent unit
    precision_tp=116, precision_fp=4,    # 120 detections checked by hand
    recall_tp=68, recall_fn=33,          # 101 real objects located independently
    reported_value=26_820,               # what the register reports
)

result = correct_unit(unit)
print(f"{result.mean:,.0f}  99% CI {result.ci(0.99)}")
print(evaluate(result, unit.reported_value).status)   # below | within | above

The verdict is decided by the interval, not the gap. A large discrepancy inside a wide interval is not a finding; a small one outside a tight interval is.

The paper's data

from bayesian_pv_census import correct_batch, load_demo

units = load_demo()                      # 93 French reporting units, kWp
correct_batch(units, prior="empirical_bayes")["corrected_mean"].sum()
# 4_033_003 kWp, the paper's national estimate

load_demo() returns the 93 continental French reporting units behind the published audit: detected rooftop capacity below 36 kWp, the manual validation counts (31,853 annotations), the size-weighted rates of specification B, and the reference values under audit from the French transmission system operator's grid-connection registry, whose redistribution RTE cleared for this release.

The full battery reproduces the paper: 25 units under-reported and 8 over-reported under specification A, an 18-unit hard core and a 7-unit negative control.

Installation counts are deliberately absent for a different reason. The audit in the paper is about capacity; shipping a count column would imply a quantity that was never audited.

What the correction assumes

The estimator is raw × P/R. Four assumptions stand behind it, and the package can only speak to two of them.

H4 binds only when the quantity is a size. If you correct a count of objects rather than a capacity, H4 drops out entirely — which is why the field is named raw and not raw_capacity.

Two of the four are therefore assumptions you carry, not results the package delivers. Reporting a credible interval without saying which of H2 and H3 you believe, and why, states less than it appears to.

Before you annotate

The interval's half-width has a closed form, so the annotation effort can be budgeted in advance rather than discovered afterwards:

from bayesian_pv_census import required_sample_size

b = required_sample_size(target_half_width=0.15, expected_precision=0.85,
                         expected_recall=0.65, level=0.99)
print(b.n_precision, b.n_recall)   # annotations needed per unit

Specifications and the hard core

A verdict that only holds under one way of computing precision and recall is not a finding. run_battery runs the audit under several and keeps what survives all of them.

from bayesian_pv_census import run_battery

battery = run_battery(units, prior="empirical_bayes")
battery.hard_core("below")       # flagged under-reported by every specification
battery.negative_control()       # flagged the other way, unanimously — the control group
battery.concordance_table()      # crosstab; empty off-diagonal corners mean no sign flips

Specification A counts annotated objects. Specification B weights them by size. B is implemented exactly as in the paper, which computes the weighted rates as point estimates and applies them as a deterministic rescaling of A's posterior — so B's interval is A's interval, shifted. That is a known limitation of the published method; the package reproduces it rather than improving on it, because reproducing the paper is the point. See the docstring of correct_unit for what a properly weighted posterior would require, and for why the min_weighted_n fallback matters more than its name suggests.

A third axis needs no specification of its own. Rescaling every raw quantity by a constant — a different surface-to-power coefficient, a different filtering threshold — multiplies the posterior and both its bounds while the reported value stays put, so it is passed generically:

run_battery(units, specs={"A": {}, "B": {"weighted": True},
                          "C_low": {"scale": 5.5 / 5.0}})

The status is monotone in that constant and flips once, so the flipping point has a closed form and no sweep is needed:

from bayesian_pv_census import conversion_threshold
conversion_threshold(result, unit.reported_value, "below")

The coefficient itself stays outside the engine. raw arrives already converted, and a conversion_coefficient argument would import photovoltaics into a core that knows nothing about it.

What forming the factor at a coarser scale costs

Correcting each unit and summing is not the same as pooling the units' rates and correcting once. 1/R is convex, so pooling always yields a smaller factor and a smaller total. The cost is exactly zero when recall is homogeneous across the pooled units, whatever the dispersion of precision, and second order in the coefficient of variation of recall otherwise:

from bayesian_pv_census import compare_aggregation_scales, pooling_penalty

compare_aggregation_scales(units, groups={"north": [...], "south": [...]}).totals
pooling_penalty(units)   # CV(R)^2 - rho * CV(P) * CV(R), and its two terms

The practical consequence is that there is no optimal grid to search for. Forming the factor at the finest level your annotation budget supports is weakly better in every case.

Scope

Out of scope by design, and unlikely to change: geospatial sampling, annotation tooling, temporal alignment between imagery and register, and detector-specific parsers. Those are format- and project-specific; this package starts once you can write down a raw total and a validation sample.

Tests

pip install -e '.[dev]' && pytest

The suite validates the mathematics on synthetic data — the closed form against the bootstrap, the monotonicity of the budgeting rule, the structural invariants of the hard core, the exact vanishing of the pooling penalty under homogeneous recall — and then checks the published French numbers against the shipped demo data, so a regression in the engine cannot pass silently.

Citation

See CITATION.cff. Licence MIT.