.. nb:name: cs6 .. nb:database: cs6 Case Study: City Air-Quality Sensor Network =========================================== This case study applies ProvSQL's continuous distributions (see :doc:`the chapter on continuous distributions `) to an air-quality sensor network, through ProvSQL Studio (see :doc:`the Studio chapter `), whose plots of densities, histograms and circuits suit random variables better than text output. The Scenario ------------ A municipal observatory runs a small network of sensors reporting a :math:`\mathit{PM}_{2.5}` concentration (fine particulate matter, in micrograms per cubic metre). Each kind of sensor has its own noise model: a normal distribution for high-end units, uniform for low-cost ones, exponential for a drift-prone unit, Erlang for a multi-pass one; a reference station reports exact values. Readings fall into three categories, *Good*, *Moderate* and *Unhealthy* (a simplified form of the US EPA breakpoints), each station is in calibration with some probability, and yesterday's readings are kept in a table of the same shape. Your tasks: * compute the probability that a reading exceeds a threshold, and its expectation beyond it; * model calibration as a mixture, aggregate readings per district, and watch the simplifier fold the result; * combine batches with ``UNION``, and compare probability methods; * find the worst reading per district, percentiles, and expected exceedances; * couple two sensors through a shared plume, and measure the coupling; * compute fleet statistics when sensors drop out; * read calibration scores as graded confidence. Setup ----- .. nb:skip .. tip:: **Prefer not to install? Use the Playground.** Instead of the manual setup below, run this Studio-driven case study in full in the ProvSQL Playground: open it `as a runnable notebook `_, or open the bare `cs6 database `_ and follow the Studio steps as you read. See the :ref:`Playground note `. .. nb:omit-begin This case study assumes a working ProvSQL installation (see :doc:`getting-provsql`) and a running ProvSQL Studio session pointed at it (see :doc:`studio`). Download :download:`setup.sql <../../casestudy6/setup.sql>` and load it into a fresh PostgreSQL database: .. code-block:: bash createdb air_quality_demo psql -d air_quality_demo -f setup.sql .. nb:omit-end .. nb:setup: ../../casestudy6/setup.sql The script creates the schema below and seeds the random-variable readings via the constructors documented in :doc:`the chapter on continuous distributions `. It is five data tables plus a provenance mapping: * ``stations(id, name, district)`` -- four monitoring stations across two districts, provenance-tracked. * ``readings(station_id, ts, pm25 random_variable)`` -- one ``pm25`` reading per station per timestamp; the ``random_variable`` carries the per-station noise model (normal, uniform, exponential, erlang, or a deterministic lifted from the reference station). * ``calibration_status(station_id, p)`` -- Bernoulli probability that each station is in calibration on the day of interest. * ``categories(name, lo, hi)`` -- three regulatory categories (*Good* / *Moderate* / *Unhealthy*) keyed by their interval bounds. * ``historical_readings(...)`` -- same shape as ``readings``, populated from yesterday's batch. * ``station_mapping(value, provenance)`` -- a provenance mapping from each reading's token to its station's name, so that :sqlfunc:`sr_formula` and PROV-XML export show station names instead of raw UUIDs. .. nb:omit-begin Connect Studio to the fixture: .. code-block:: bash provsql-studio --dsn postgresql:///air_quality_demo and open `http://127.0.0.1:8000/ `_ in a browser. .. nb:omit-end The schema panel lists the fixture's six relations: the four provenance-tracked tables (``stations``, ``calibration_status``, ``readings``, ``historical_readings``) carry the purple :sc:`prov` pill, ``categories`` is plain, and ``station_mapping`` is tagged :sc:`mapping`. The ``pm25`` column on ``readings`` and ``historical_readings`` is flagged with a terracotta :sc:`rv` pill: comparison and arithmetic operators on this column build provenance circuits, so a condition like ``pm25 > 35`` yields an uncertain event, not a Boolean. .. figure:: /_static/casestudy6/schema-panel.png :alt: Studio schema panel listing readings, historical_readings, calibration_status, stations (all PROV-tagged), categories (no provenance), and the station_mapping table, with the pm25 column on readings and historical_readings flagged with a small RV pill. The schema panel opened from the top nav. The four provenance-tracked tables carry the purple :sc:`prov` pill; ``readings`` and ``historical_readings`` list ``pm25`` with a terracotta :sc:`rv` pill marking it as a ``random_variable`` column. .. _cs6-step-1: Step 1: Inspect a Noisy Reading -------------------------------- In the Studio query box: .. code-block:: postgresql SELECT id, ts, pm25 FROM readings WHERE station_id = 's1' ORDER BY ts The result table renders ``pm25`` as a clickable ``random_variable`` cell carrying the underlying gate UUID. Click into a row's ``pm25``: Studio renders its circuit, a single ``gate_rv`` leaf with the distribution-kind initial in the circle (N for a Normal, U for Uniform, Exp for Exponential, Erl for Erlang). Pick :guilabel:`Distribution profile` from the :guilabel:`Distribution` group of the eval strip and click :guilabel:`Run`: the panel returns headline stats: support, mean :math:`\mu`, standard deviation :math:`\sigma`, and entropy :math:`H`, as well as an inline histogram with a PDF/CDF toggle. The histogram is backed server-side by :sqlfunc:`rv_histogram`; pinning ``provsql.monte_carlo_seed`` in the Config panel (under :guilabel:`Evaluation`) makes the shape reproducible across re-runs. .. figure:: /_static/casestudy6/gate-rv-distribution-profile.png :alt: Studio Circuit mode showing the gate_rv N(28,2) leaf at the top of the canvas and the Distribution profile eval-strip panel below, with mu=28, sigma=2, support (-infinity, +infinity) and an inline histogram overlaid by a terracotta analytical PDF curve. The ``gate_rv`` leaf for ``pm25`` on row 1 is a ``N(28, 2)`` circle; the eval-strip :guilabel:`Distribution profile` panel shows the :math:`\mu`, :math:`\sigma`, :math:`H`, support, and an inline histogram. When the gate has a closed-form family (here a bare Normal), Studio overlays the analytical PDF on top of the bars. .. _cs6-step-2: Step 2: A First Probabilistic Threshold ---------------------------------------- The :math:`\mathit{PM}_{2.5}` *Unhealthy* category begins above 35.1. Find the rows whose reading might cross it: .. code-block:: postgresql SELECT id, station_id, ts FROM readings WHERE pm25 > 35 Because ``pm25`` is a random variable, the comparison is not a yes/no test: it stands for the *event* "this reading exceeds 35", and ProvSQL attaches that event to each row's provenance. Click into a result row's auto-added ``provsql`` cell. Circuit mode shows the Boolean wrapper (a ``gate_times`` over the row's input token and the ``gate_cmp``); the cmp's child link reaches into the ``gate_rv`` from :ref:`Step 1 `. The eval strip's :guilabel:`Marginal probability` entry exposes the full method catalogue (see :doc:`the chapter on probabilities `). Pick ``monte-carlo`` and set ``n = 10000``; the panel returns the probability with a Hoeffding confidence band. Pin ``provsql.monte_carlo_seed = 42`` in the Config panel and re-run: the result is now identical across runs. Set the seed back to ``-1`` and re-run to see the band shift between runs. .. figure:: /_static/casestudy6/cmp-circuit-mc-eval.png :alt: Studio Circuit mode showing the gate_times wrapper above the iota input gate (left) and the > gate_cmp (right); the cmp's children are the N(28,2) gate_rv and the constant 35. The provenance of one row from ``WHERE pm25 > 35``: a ``gate_times`` (``⊗``) wraps the row's input token ``ι`` and a ``gate_cmp`` ``>`` whose children are the ``N(28, 2)`` leaf and the constant ``35``. The eval strip below switches to ``probability_evaluate`` and exposes the method picker. Step 3: The Simplifier in Action --------------------------------- Each comparison is recorded as a ``gate_cmp``. A *simplifier* then folds comparators whose answer can be decided from the operands' support alone: for example, ``U(10, 22) > 35`` is universally false because the uniform's upper bound is below the threshold. The fold is controlled by ``provsql.simplify_on_load`` (default on), which the Config panel exposes under :guilabel:`Evaluation`. Click into row 2's auto-added ``provsql`` cell from the :ref:`Step 2 ` result (station ``s2``, ``pm25 ~ U(10, 22)``). With ``provsql.simplify_on_load`` on, the canvas shows a single ``𝟘`` (zero) gate: the simplifier resolved the comparator to a constant-false leaf and dropped the whole subtree. Toggle the GUC off in the Config panel and click the cell again: the canvas now shows the raw construction shape, a ``gate_times`` (``⊗``) over the row's input token ``ι`` and a ``gate_cmp`` (``>``) whose children are the ``U(10, 22)`` leaf and the constant ``35``. Both views are semantically identical. .. figure:: /_static/casestudy6/simplify-before-after.png :alt: Side-by-side: on the left, the raw circuit for row 2's provenance with provsql.simplify_on_load off, showing a gate_times over an input token and a > gate_cmp whose children are U(10,22) and 35; on the right, the same circuit with the GUC on collapsed to a single zero gate. Row 2's provenance for ``pm25 > 35`` with ``provsql.simplify_on_load`` toggled off (left) vs on (right). The simplifier recognised that the upper bound of ``U(10, 22)`` is below the threshold, so the comparator is universally false and the whole subtree collapses to the additive identity ``𝟘``. Step 4: Calibration via Mixtures --------------------------------- Each station has a probability of being mis-calibrated; a mis-calibrated unit over-reports by 20% (the *reading* it records is ``1.2`` times the true value). The corrected estimate of the true reading is therefore ``pm25`` with probability ``p`` (the station is in spec) and ``pm25 / 1.2`` with probability ``1 - p`` (the report needs to be scaled back). Express this as a Bernoulli mixture: .. code-block:: postgresql SELECT r.id, r.station_id, mixture(cs.p, r.pm25, r.pm25 / 1.2) AS pm25_calibrated FROM readings r JOIN calibration_status cs USING (station_id) WHERE r.station_id = 's1' Click into a result row's ``pm25_calibrated`` cell. Circuit mode renders the ``gate_mixture`` as a ``Mix`` node with three labelled outgoing edges (``p`` / ``x`` / ``y``) matching the SQL constructor's argument order: ``p`` points to the Bernoulli mixing probability, ``x`` to the in-spec arm, and ``y`` to the correction arm. The same node-inspector panel exposes ``Distribution profile`` on the mixture root. Because station ``s1`` is in spec 95% of the time, the histogram is dominated by the ``N(28, 2)`` arm and the out-of-spec ``N(23.33, 1.667)`` contributes only a small left shoulder rather than a visually distinct second mode; the panel headline reflects this with a mixture mean slightly below 28. To see clear bimodality, re-run the query with a larger calibration error, e.g., replace ``r.pm25 / 1.2`` with ``r.pm25 / 2.0`` so the out-of-spec arm folds to ``N(14, 1)``, well separated from the in-spec ``N(28, 2)``; the two peaks then show up distinctly on the histogram even at the 95%/5% weighting. .. figure:: /_static/casestudy6/mixture-node.png :width: 696px :alt: Mix node with three labelled outgoing edges (p, x, y); the p child is a 95% Bernoulli, the x child is the N(28,2) reading, and the y child is the back-scaled N(23.33,1.667) folded by the simplifier. The ``gate_mixture`` for the calibrated reading. The ``95%`` child is the Bernoulli probability that station ``s1`` is in spec; the ``x`` arm is the raw reading ``N(28, 2)``; the ``y`` arm is the back-scaled estimate ``pm25 / 1.2``, which the simplifier folded through the Normal affine-shift rule into a single ``gate_rv`` ``N(23.33, 1.667)``. Circle labels show four significant figures; the inspector pinned to either child surfaces the full-precision parameters. .. _cs6-step-5: Step 5: Aggregation Over Random Variables ------------------------------------------ Compute average :math:`\mathit{PM}_{2.5}` per district: .. code-block:: postgresql SELECT s.district, avg(r.pm25) AS avg_pm25, sum(r.pm25) AS total_pm25 FROM readings r JOIN stations s ON s.id = r.station_id GROUP BY s.district Click into a row's ``avg_pm25`` cell. Circuit mode shows the :sqlfunc:`avg` circuit: a ``gate_arith(DIV, num, denom)`` over two ``gate_arith(PLUS, …)`` subtrees, each child a per-row ``gate_mixture``. The right child of the outer division is the count of *included* rows under their per-row provenance: rows whose provenance is false contribute the additive identity to both numerator and denominator. Run :guilabel:`Distribution profile` on the root: the panel shows the per-district average as a tight distribution centred at the inclusion-weighted mean. .. figure:: /_static/casestudy6/sum-of-mixtures.png :width: 596px :alt: DAG with a single root division node, two PLUS subtrees under it, eight Mix nodes at the next level, and the gate_rv leaves (N(28,2), N(40,4), U(12,24), U(10,22)) plus iota input gates at the bottom. The ``avg(pm25)`` cell for the *centre* district lowers to ``gate_arith(DIV, gate_arith(PLUS, mixtures), gate_arith(PLUS, one-mixtures))``. The eight ``Mix`` nodes are the district's four readings (two stations × two timestamps), each appearing once under the numerator's sum and once under the denominator's count; the ``gate_rv`` leaves at the bottom are the per-reading distributions; the ``ι`` leaves anchor each row's provenance. .. _cs6-step-6: Step 6: Conditional Inference ------------------------------ Re-open the filtered query from :ref:`Step 2 `: .. code-block:: postgresql SELECT id, station_id, ts, pm25 FROM readings WHERE pm25 > 35 AND station_id = 's1' Click a result row's ``pm25`` cell and pick :guilabel:`Distribution profile` in the eval strip. The strip then shows a :guilabel:`Condition on` text input, preset to the row's provenance UUID, with the :guilabel:`Conditioned by:` badge underneath it active. Run: the histogram now shows the *truncated* shape, restricted to the tail above ``35``. Pick :guilabel:`Moment` with ``k = 1`` and ``raw``: the panel returns the closed-form Mills-ratio mean of the `truncated normal `_, exactly :math:`\mu + \sigma \cdot \frac{\phi(\alpha)}{1 - \Phi(\alpha)}` with :math:`\alpha = (35 - \mu)/\sigma`. Click the active badge to clear the conditioning; the panel reverts to the unconditional mean :math:`\mu`. Click the muted badge to restore the row provenance. The closed-form truncation table covers every family implementing truncated moments -- Normal (Mills ratio), Uniform (intersected support), Exponential (memorylessness on a lower bound or finite-interval truncation), Log-normal, Weibull, Pareto (tail self-similarity), and Beta. For other shapes, the conditional moment is estimated by rejection sampling at budget ``provsql.rv_mc_samples``. .. figure:: /_static/casestudy6/condition-on-active.png :width: 696px :alt: The Distribution profile panel showing supp [35, +infinity], mu approximately 40.82, sigma approximately 3.35; the Condition on input is populated with the row's provenance UUID and the Conditioned by badge is active. The conditional distribution profile for row 5 (``pm25 ∼ N(40, 4)``) under the event ``pm25 > 35``. Studio auto-presets the :guilabel:`Condition on` input with the row's provenance UUID and activates the :guilabel:`Conditioned by:` badge; the panel's header reflects the truncated support ``[35, +∞]`` and the Mills-ratio mean ``μ ≈ 40.82``, ``σ ≈ 3.35`` (closed form on the truncated normal). Step 7: Diagnostic Sampling ---------------------------- For raw inspection or downstream analytics, draw samples from the conditional distribution. With the same row pinned and the :guilabel:`Conditioned by:` badge active, pick :guilabel:`Sample` from the :guilabel:`Distribution` group; set ``n = 200`` and run. The result panel shows a six-value inline preview with a "show full list" expander; clicking it dumps all 200 samples. Both rows of the query are bare Normals conditioned on an interval, so all 200 samples come straight from the truncated distribution; for other shapes, :sqlfunc:`rv_sample` falls back to rejection sampling and may return fewer (see :doc:`continuous-distributions`). Step 8: Combining Batches via UNION ------------------------------------ Both batches share the same id space (rows ``1`` through ``8``, one per ``(station, timestamp)`` slot), so a ``UNION`` (without ``ALL``) over ``(station_id, id)`` deduplicates a slot to a single result row whose provenance combines today's reading and yesterday's reading via the semiring addition. With ``WHERE pm25 > 35`` lifted on each branch, each contributing row carries a ``gate_cmp(pm25 > 35)`` of its own: .. code-block:: postgresql (SELECT station_id, id FROM readings WHERE pm25 > 35) UNION (SELECT station_id, id FROM historical_readings WHERE pm25 > 35) ORDER BY station_id, id Pick a result row's ``provsql`` cell: Circuit mode shows a ``gate_plus`` (``⊕``) over the two contributing inputs (today's row and yesterday's row), each carrying its own ``gate_cmp(pm25 > 35)`` from the lifted ``WHERE``. Run :guilabel:`Marginal probability` from the eval strip: it gives the probability that *at least one* of the two days produced an Unhealthy reading for that slot. We deliberately keep the ``random_variable`` ``pm25`` column out of the ``SELECT``: there is no duplicate-elimination semantics for ``random_variable``, so a ``UNION`` over an RV column would have no well-defined meaning. Step 9: Filtering Grouped Random Variables by Expected Value ------------------------------------------------------------- Filter the per-district aggregates from :ref:`Step 5 ` by their expected average. Because :sqlfunc:`avg` over a ``random_variable`` column returns a ``random_variable`` (not an ``agg_token``), and :sqlfunc:`expected` collapses it to a plain ``double``, the ``HAVING`` condition is an ordinary comparison that PostgreSQL evaluates directly; each surviving group's provenance still gets a ``delta(gate_agg)`` wrapper: .. code-block:: postgresql SELECT s.district, avg(r.pm25) AS avg_pm25 FROM readings r JOIN stations s ON s.id = r.station_id GROUP BY s.district HAVING expected(avg(r.pm25)) > 20 :sqlfunc:`expected` collapses the inner :sqlfunc:`avg` (the division circuit of :ref:`Step 5 `) to its mean. It does so exactly here: no probability has been written on these rows' tokens, so every row is certainly present, the count in the denominator is a constant, and the mean of the average is the average of the means. ``> 20`` then compares two ``double`` values, so a group survives exactly when its expected average exceeds the threshold. For the case-study fixture both districts pass (centre at 25.5, east at 21.625); clicking either result row's ``provsql`` cell shows the ``delta(gate_agg)`` shape, identical to the no-HAVING aggregate from :ref:`Step 5 ` but filtered to the surviving groups. Step 10: Independent vs Monte Carlo ------------------------------------ For threshold queries whose contributing rows have structurally independent provenance, the ``'independent'`` probability method (see :doc:`the chapter on probabilities `) is *exact* and far cheaper than Monte Carlo. Re-run the :ref:`Step 2 ` query, click a result row's ``provsql`` cell, and run :guilabel:`Marginal probability` from the eval strip three times, with the methods ``independent``, ``tree-decomposition``, and ``monte-carlo`` at ``n = 10000``: the two exact methods return the same value to full precision, while ``monte-carlo`` returns a Hoeffding-bounded estimate that tightens as ``n`` grows. From SQL, the same methods are arguments of :sqlfunc:`probability_evaluate`. Step 11: The Worst Reading in a District ----------------------------------------- Which district had the worst air this morning? The question asks for the distribution of the *maximum* over each district's readings -- an order statistic, not a sum. The ``max`` aggregate over a ``random_variable`` column builds exactly that: .. code-block:: postgresql SELECT s.district, max(r.pm25) AS worst, expected(max(r.pm25)) AS worst_mean FROM readings r JOIN stations s ON s.id = r.station_id GROUP BY s.district ORDER BY s.district The *centre* district's worst reading averages ≈ 40.0 -- its ``N(40, 4)`` reading dominates the maximum -- and *east* comes out at ≈ 39.6 (a Monte Carlo estimate, see below, so it varies by a few tenths between runs and with the seed), driven by the heavy right tails of ``Exp(0.04)`` and ``erlang(3, 0.1)`` even though both means sit near 25--30. Note how the extremum tells a different story from the averages of :ref:`Step 5 ` (≈ 25.5 and ≈ 21.6): the east district looks fine on average and just as alarming at the extreme. ``worst_mean`` itself is a plain number -- ``expected`` has already collapsed the distribution to its mean -- which is why the query also projects the unreduced aggregate: the ``worst`` cell is a clickable ``random_variable`` token whose circuit is a ``gate_arith`` ``MAX`` root over the per-row mixtures; a row absent in a world contributes ``-∞``, the order-statistic identity, so it never perturbs the maximum. Since every row is certainly present, :sqlfunc:`expected` reads the maximum of the readings themselves: over continuous distributions, as in the centre district, it is computed exactly by a one-dimensional integral; the east district's constant reference readings fall outside that formula, and its expectation is estimated by Monte Carlo. The same-row form is ``greatest`` / ``least``: the SQL keywords accept ``random_variable`` arguments, which ProvSQL rewrites into order statistics (:ref:`Step 15 ` uses ``least`` outside any provenance-tracked relation): .. code-block:: postgresql SELECT greatest(a.pm25, b.pm25) AS worse, expected(greatest(a.pm25, b.pm25)) AS worse_of_two FROM readings a, readings b WHERE a.id = 1 AND b.id = 2 Clicking the ``worse`` cell shows the same ``gate_arith`` ``MAX`` root, but here the children *are* independent bare leaves (``N(28, 2)`` and ``U(10, 22)``), so the mean is computed analytically by the layer-cake integral -- ≈ 28.00003, barely above the normal's mean, since the uniform almost never wins. Step 12: Percentiles and Exceedance Margins -------------------------------------------- :sqlfunc:`quantile` reads the inverse CDF off any reading -- the natural regulatory question is "what level is only exceeded 5% of the time?": .. code-block:: postgresql SELECT id, station_id, pm25, quantile(pm25, 0.95) AS p95 FROM readings WHERE id IN (1, 3, 5, 7) ORDER BY id Each family answers through its own inverse CDF: the normals give ``μ + 1.645·σ`` (≈ 31.29 for ``N(28, 2)``, ≈ 46.58 for ``N(40, 4)``), the exponential gives ``−ln(0.05)/λ ≈ 74.89``, and the Erlang -- which has no elementary inverse -- is bisected on its closed-form CDF to ≈ 62.96. No sampling is involved in any of these. The same readout is available interactively: with a ``pm25`` cell pinned, pick :guilabel:`Quantile` from the :guilabel:`Distribution` group of the eval strip, set the fraction ``p``, and run. Quantiles compose with conditioning: the median of a reading *given* that it crossed the Unhealthy line is the ``p = 0.5`` quantile of the truncated distribution, again in closed form: .. code-block:: postgresql SELECT id, pm25, quantile(pm25, 0.5, provenance()) AS median_given_high FROM readings WHERE pm25 > 35 AND station_id = 's1' ORDER BY id Row 1 (``N(28, 2)``, for which the event is a far tail) gives ≈ 35.36, just above the threshold; row 5 (``N(40, 4)``) gives ≈ 40.53. In the eval strip, the same conditioning rides on the :guilabel:`Conditioned by:` badge from :ref:`Step 6 `: with it active, :guilabel:`Quantile` returns the truncated-distribution quantile. Step 13: Expected Excess via CASE ---------------------------------- *How far* above the threshold does a station land, on average? The excess ``max(pm25 − 35, 0)`` is a piecewise transform, written as a searched ``CASE`` (see :doc:`the continuous distributions chapter `): .. code-block:: postgresql SELECT id, expected(CASE WHEN pm25 > 35 THEN pm25 - 35 ELSE as_random(0) END) AS expected_excess FROM readings WHERE station_id = 's1' ORDER BY id The guard and the branches read the same ``pm25``, and a piecewise function of one random variable is among the ``CASE`` shapes whose moments come back in closed form: the expectation integrates each branch over its guard's region, so the correlation between "the reading exceeded 35" and "by how much" is preserved exactly, with no sampling -- the answers hold even under ``SET provsql.rv_mc_samples = 0``. Row 1 returns ≈ 0.0001 -- ``N(28, 2)`` essentially never crosses -- while row 5 returns the closed-form partial expectation :math:`\sigma\,\phi(\alpha) + (\mu - 35)\,\Phi(-\alpha)` with :math:`\alpha = (35 - \mu)/\sigma`, ≈ 5.2023 for ``N(40, 4)``, to full precision. Note the explicit ``as_random(0)`` (or equivalently ``0::random_variable``) on the ``ELSE`` branch: ``CASE`` type resolution needs the branches in a single type category, so a bare numeric literal will not lift on its own. .. _cs6-step-14: Step 14: Correlated Sensors ---------------------------- The two centre-district stations sit a few hundred metres apart; on a windy day both are shifted by the same dust plume. Model the plume once and add it to both readings -- the shared leaf is what makes them correlated -- then measure it: .. code-block:: postgresql WITH plume AS (SELECT normal(0, 3) AS dust) SELECT covariance(a.pm25 + p.dust, b.pm25 + p.dust) AS cov_shared, correlation(a.pm25 + p.dust, b.pm25 + p.dust) AS corr_shared, covariance(a.pm25, b.pm25) AS cov_indep, stddev(a.pm25) AS sd_a FROM readings a, readings b, plume p WHERE a.id = 1 AND b.id = 2 The two raw readings are structurally independent -- their base-RV footprints are disjoint -- so ``cov_indep`` is *exactly* ``0`` (no sampling), and ``sd_a`` is exactly ``2``. The plume-shifted pair shares the ``dust`` leaf: in theory ``Cov = Var(dust) = 9`` and ``ρ = 9 / √((4+9)·(12+9)) ≈ 0.545``; the readouts land close (≈ 9.06 and ≈ 0.544 at seed 42). A shared-leaf cross-moment has no closed form, so both statistics are estimated by Monte Carlo from coupled draws, the plume being drawn once per iteration and shared by both readings, which keeps the estimates tight around the theoretical values. .. _cs6-step-15: Step 15: Drift, Lifetime, and Spikes: More Families ---------------------------------------------------- The vendor spec sheets bring in three more families, all available as constructors (see :doc:`the chapter on continuous distributions `). The multi-pass unit's zero drift accumulates as a Gamma; its per-batch drift check is a chi-squared test statistic: .. code-block:: postgresql SELECT expected(gamma(2.5, 0.5)) AS drift_mean, variance(gamma(2.5, 0.5)) AS drift_var, quantile(chi_squared(4), 0.95) AS chi2_crit ``gamma(k = 2.5, λ = 0.5)`` has mean ``k/λ = 5`` and variance ``k/λ² = 10``, both exact; the χ²₄ critical value comes out ≈ 9.4877, bisected on the regularised-incomplete-gamma CDF. Sensor lifetimes are Weibull. A station carrying two redundant units fails over to the second when the first dies; the time to *first* failure is the minimum of the two lifetimes: .. code-block:: postgresql SELECT expected(least(weibull(1.5, 4.0), weibull(1.5, 4.0))) AS first_failure, expected(weibull(1.5, 4.0)) AS single_lifetime A single ``weibull(k = 1.5, λ = 4)`` unit lasts ``λ·Γ(1 + 1/k) ≈ 3.61`` years on average; the minimum of two i.i.d. Weibulls is again Weibull (min-stability), so the first failure at ``≈ 2.27`` years is computed in closed form. Construction-season episodic spikes are heavy-tailed -- a Pareto over a 30 μg/m³ floor: .. code-block:: postgresql SELECT expected(pareto(30, 2.5)) AS spike_mean, quantile(pareto(30, 2.5), 0.95) AS spike_p95, expected(pareto(30, 1.0)) AS undefined_mean The ``α = 2.5`` spike averages ``α·xₘ/(α−1) = 50`` with a 95th percentile near 99.4 -- the heavy tail at work. The last column shows moment honesty: ``pareto(30, 1)`` has *no* finite mean, and ProvSQL reports ``Infinity`` rather than estimating a divergent integral. Finally, does a spike beat the noisy high reading? A Pareto vs Normal comparison has no registered closed form, so ProvSQL integrates the pair by quadrature -- exact-grade, no sampling: .. code-block:: postgresql SELECT probability(x > y) AS p_spike_beats_normal FROM (SELECT pareto(30, 2.5) AS x, normal(45, 5) AS y) t which returns ≈ 0.384: even though the spike's *mean* (50) is above the normal's (45), the spike spends most of its mass near its 30 floor and only wins through its tail. Step 16: Fleet Statistics Under Sensor Dropout ----------------------------------------------- Every reading so far contributed unconditionally. Today's maintenance log says otherwise: a station out of calibration silently drops out of the feed. Pin each reading's provenance probability to its station's calibration probability, so a row is *present* only when its station is in calibration: .. code-block:: postgresql SELECT set_prob(provenance(), (SELECT p FROM calibration_status cs WHERE cs.station_id = r.station_id)) FROM readings r The calibration probability is fetched through a scalar subquery so that ``provenance()`` stays the ``readings`` row's own input token: joining the two tracked tables instead would make it the *joint* row provenance (a ``times`` gate), which ``set_prob`` cannot pin. ProvSQL warns that the subquery's data is read outside provenance tracking ("treated as certain") -- exactly the intent here: ``p`` parameterises the presence model, it is not uncertain data itself. The SQL-standard statistic aggregates (:sqlfunc:`stddev_pop` / :sqlfunc:`stddev_samp`, :sqlfunc:`covar_pop` / :sqlfunc:`covar_samp` / :sqlfunc:`corr`, and the ordered-set :sqlfunc:`percentile_cont`) are lifted to ``random_variable`` rows with exactly these semantics: a row absent in a world drops out of every sum, the count, and the percentile's member set. Start where the answer is checkable by hand -- the reference station is always in calibration (``p = 1``) and its readings are lab-grade constants ``15.0`` and ``16.5``: .. code-block:: postgresql SELECT expected(stddev_pop(pm25)) AS sd_ref, expected(corr(pm25, pm25 * 2)) AS corr_ref FROM readings WHERE station_id = 's4' Both statistics are Dirac -- every Monte-Carlo draw is the same world and the same constants -- so the answers are exact: ``sd_ref = 0.75`` (the population stddev of ``{15, 16.5}``) and ``corr_ref = 1`` (a reading is perfectly correlated with its own rescaling, drawn once and observed by both argument positions). Now the probabilistic fleet: the *median* reading per district, where dropout changes the member set itself -- when Riverside Park is out of calibration (30% of days), both of its rows leave the group: .. code-block:: postgresql SELECT s.district, expected(percentile_cont(0.5) WITHIN GROUP (ORDER BY r.pm25)) AS median_pm25 FROM readings r JOIN stations s ON s.id = r.station_id GROUP BY s.district In each Monte-Carlo draw, only the values whose station is in calibration are kept, sorted, and interpolated (pin the result node and look for the ``p50`` circle in the canvas). Under dropout the centre district's expected median lands around 26 -- pulled between the two Normal City-Centre readings and the two lower Riverside uniforms that are only present 70% of the time. Step 17: A Maintenance-Triage Headline via CASE over Aggregates ---------------------------------------------------------------- Which district needs the maintenance crew first? The desk wants a single *headline* number per district: if even the weakest station's calibration confidence clears 0.8, report the district's total confidence weight; otherwise flag the weakest station (the crew's target). That is a ``CASE`` whose guards **and** branches are aggregates (see :doc:`the aggregation chapter `). Its expectation is exact for ``sum`` / ``count`` / ``min`` / ``max`` guards and branches; an ``avg`` branch inside a ``CASE`` is exact only for an unconditional expectation and is evaluated by Monte Carlo otherwise. First make the maintenance log itself uncertain (each record is confirmed with a probability) and the station registry certain: .. code-block:: postgresql SELECT set_prob(provenance(), 1.0) FROM stations; SELECT set_prob(provenance(), CASE station_id WHEN 's1' THEN 1.0 WHEN 's2' THEN 0.8 WHEN 's3' THEN 0.7 ELSE 1.0 END) FROM calibration_status; One record per district is kept certain, so each district's group exists in every world (an all-absent group would produce no headline row at all). .. code-block:: postgresql SELECT s.district, CASE WHEN min(cs.p) > 0.8 THEN sum(cs.p) ELSE min(cs.p) END AS confidence_headline FROM calibration_status cs JOIN stations s ON s.id = cs.station_id GROUP BY s.district The result cells display as ``0.7 (*)`` and ``0.6 (*)`` -- the ``CASE`` evaluated on the actual data (with every record confirmed, both districts have a weak link below the bar, so each headlines its worst station), with the ``(*)`` marking a probabilistic aggregate whose value is world-dependent. The headline's distribution is evaluated *exactly*, with no Monte Carlo, even under ``SET provsql.rv_mc_samples = 0``: .. code-block:: postgresql SELECT s.district, expected(CASE WHEN min(cs.p) > 0.8 THEN sum(cs.p) ELSE min(cs.p) END) AS expected_headline FROM calibration_status cs JOIN stations s ON s.id = cs.station_id GROUP BY s.district For the centre district there are exactly two worlds: with both records confirmed (probability ``0.8``) the weakest link ``0.70`` misses the bar and is the headline; with only City Centre's record confirmed (``0.2``) the sole confidence ``0.95`` clears it and the total ``0.95`` is reported. So :math:`E = 0.8 \cdot 0.70 + 0.2 \cdot 0.95 = 0.75` exactly -- and that is the value the query returns. The east district works out to :math:`0.7 \cdot 0.60 + 0.3 \cdot 1.00 = 0.72`. Step 18: Shared Information in Nats ------------------------------------ :ref:`Step 14 ` measured the plume coupling with :sqlfunc:`correlation`; :sqlfunc:`mutual_information` is its information-theoretic counterpart -- symmetric, in nats, and zero *exactly* when the two variables share no randomness: .. code-block:: postgresql WITH plume AS (SELECT normal(0, 3) AS dust) SELECT mutual_information(a.pm25, b.pm25) AS mi_indep, mutual_information(a.pm25 + p.dust, b.pm25 + p.dust) AS mi_shared FROM readings a, readings b, plume p WHERE a.id = 1 AND b.id = 2 ``mi_indep`` is **exactly 0** with no sampling at all: the two raw readings share no random leaf. The plume-shifted pair shares the ``dust`` leaf, and the readout estimates their mutual information by a 2-D histogram over *coupled* joint draws (each iteration draws the shared plume once and both sums observe it) -- about ``0.21`` nats at seed 42. For a jointly-Gaussian pair with the :ref:`Step 14 ` correlation :math:`\rho \approx 0.545` the closed form would give :math:`-\tfrac12\ln(1-\rho^2) \approx 0.176`; the estimate lands in that neighbourhood, a little above: the pair is not jointly Gaussian (Riverside's noise is uniform), and the plug-in histogram estimator carries a small upward bias. Step 19: An Offline Sensor Reports NULL ---------------------------------------- Sensors go offline. When Riverside Park misses its 10:00 sample, the natural encoding is a row whose ``pm25`` is **SQL NULL** -- there is no reading, not even an uncertain one: .. code-block:: postgresql INSERT INTO readings (id, station_id, ts, pm25) VALUES (9, 's2', '2026-05-12 10:00', NULL); Three things follow, each matching plain SQL: - **Threshold events are unknown, never true.** ``WHERE pm25 > 35`` over the offline row is *unknown* under SQL's three-valued logic in every possible world, so the lifted comparison annotates the row with the circuit's zero -- :sqlfunc:`probability_evaluate` returns exactly 0, with no sampling. (A NULL reading is not "certainly exceeding", and not "possibly exceeding" either: it is no evidence at all.) - **Null tests are ordinary and deterministic.** ``WHERE pm25 IS NULL`` finds the offline row like any other predicate; PostgreSQL evaluates it and provenance simply carries the row's token through. - **Aggregates skip the offline row.** ``expected(avg(pm25))`` over Riverside Park is the same with or without row 9: the NULL reading contributes to neither the sum nor the count. Delete the placeholder row to return to the running dataset: .. code-block:: postgresql DELETE FROM readings WHERE id = 9; The general rules -- which predicates treat NULLs as unknown and what that means for the circuits -- are in :doc:`the NULL semantics chapter `. Step 20: Calibration as Graded Confidence ----------------------------------------- The calibration scores of ``calibration_status`` can also be read as *degrees of confidence*, not probabilities. A district's alert is corroborated when two of its stations agree; how confident is that corroboration? A fuzzy semiring combines the degrees, with ``max`` over alternatives and a t-norm for a conjunction. Turn the scores into a provenance mapping, and list the districts with a corroborating pair of stations: .. code-block:: postgresql DROP TABLE IF EXISTS calibration_degree; SELECT create_provenance_mapping('calibration_degree', 'calibration_status', 'p'); SELECT DISTINCT s1.district FROM calibration_status c1 JOIN stations s1 ON s1.id = c1.station_id JOIN calibration_status c2 ON c2.station_id > c1.station_id JOIN stations s2 ON s2.id = c2.station_id AND s2.district = s1.district ORDER BY s1.district; Click the ``provsql`` cell of the centre's row, pick :guilabel:`Łukasiewicz (numeric fuzzy)` from the :guilabel:`Numeric` group of the eval strip, select the ``calibration_degree`` mapping, and run; then do the same with :guilabel:`Viterbi (max-times)`. The Łukasiewicz semiring conjoins with the t-norm :math:`\max(a + b - 1, 0)`: the centre, corroborated by City Centre (0.95) and Riverside Park (0.70), gets 0.65. The two stations' shortfalls from full confidence, 0.05 and 0.30, add up. Viterbi multiplies instead, 0.665. With the fully trusted Suburban Reference (1.00), the east district keeps the Industrial Estate's 0.60 under both: a certain station costs nothing. The stations' own tokens are not in the mapping and count as fully true. In SQL, the same values are those of :sqlfunc:`sr_lukasiewicz` and :sqlfunc:`sr_viterbi`, with the mapping's name as second argument. See :doc:`the chapter on continuous distributions ` for the full reference and :doc:`the Studio chapter ` for the Studio reference.