chromy_core/
estimate.rs

1//! Estimation that matches the design.
2//!
3//! Under minimum replacement a unit may be selected more than once, and
4//! the estimator that fits is the multiplicity (Hansen–Hurwitz form)
5//! estimator: each unit's value divided by its *effective* expected hits —
6//! the ones the selection actually used — and weighted by the hits it
7//! received. For the variance there is no unbiased answer under a
8//! systematic-like design; the successive-difference estimator is the
9//! working convention, computed over the units whose selection was
10//! actually random, in list order, so that it reads the local variation
11//! the ordering left behind.
12
13use crate::error::Error;
14
15/// The totals and their successive-difference variance.
16#[derive(Debug, Clone)]
17pub struct Estimate {
18    /// The estimated total of each study variable, in column order.
19    pub total: Vec<f64>,
20    /// How many sampled units had a randomly decided hit — the sample the
21    /// variance is computed over.
22    pub random_draws: usize,
23    /// The successive-difference variance of each total: zero when no hit
24    /// was random, `NaN` when exactly one was (one observation estimates
25    /// no variance).
26    pub sd_variance: Vec<f64>,
27}
28
29/// The multiplicity estimator with successive-difference variance.
30///
31/// `y` is the matrix of study variables for the *frame*, column-major with
32/// `n_rows` units. The remaining slices describe the sample, one entry per
33/// sampled unit: `units` the 0-based frame indices, `positions` the
34/// 0-based positions in the selection list (the order the variance reads),
35/// `expected_hits` the effective expected hits, `hits` the realised hits.
36pub fn multiplicity_estimate(
37    y: &[f64],
38    n_rows: usize,
39    units: &[usize],
40    positions: &[usize],
41    expected_hits: &[f64],
42    hits: &[u64],
43) -> Result<Estimate, Error> {
44    if n_rows == 0 || y.len() % n_rows != 0 {
45        return Err(Error::Input(format!(
46            "a matrix of {} values cannot have {} rows",
47            y.len(),
48            n_rows
49        )));
50    }
51    let n_cols = y.len() / n_rows;
52    let m = units.len();
53    if positions.len() != m || expected_hits.len() != m || hits.len() != m {
54        return Err(Error::Input(
55            "units, positions, expected hits and hits must have one entry per \
56             sampled unit"
57                .into(),
58        ));
59    }
60    for &v in y {
61        if !v.is_finite() {
62            return Err(Error::Input(format!(
63                "study variables must be finite; found {v}"
64            )));
65        }
66    }
67    for &u in units {
68        if u >= n_rows {
69            return Err(Error::Input(format!(
70                "sampled unit {u} is beyond the {n_rows} rows of the study matrix"
71            )));
72        }
73    }
74    for (s, &e) in expected_hits.iter().enumerate() {
75        if !e.is_finite() || e <= 0.0 {
76            return Err(Error::Input(format!(
77                "expected hits must be positive on the sample; entry {s} is {e}"
78            )));
79        }
80    }
81
82    // The sample in list order: the successive differences are to read
83    // neighbours of the ordered frame, not of whatever order the caller
84    // stored the sample in.
85    let mut by_list: Vec<usize> = (0..m).collect();
86    by_list.sort_by_key(|&s| positions[s]);
87
88    let mut total = vec![0.0_f64; n_cols];
89    for c in 0..n_cols {
90        for &s in &by_list {
91            total[c] += y[c * n_rows + units[s]] / expected_hits[s] * hits[s] as f64;
92        }
93    }
94
95    // The units whose hit was random: realised hits one above the floor of
96    // the expected hits. Fixed hits are certainties and contribute no
97    // variance.
98    let random: Vec<usize> = by_list
99        .iter()
100        .copied()
101        .filter(|&s| hits[s] as i64 - expected_hits[s].floor() as i64 == 1)
102        .collect();
103    let draws = random.len();
104
105    let sd_variance: Vec<f64> = (0..n_cols)
106        .map(|c| {
107            if draws == 0 {
108                0.0
109            } else if draws == 1 {
110                f64::NAN
111            } else {
112                let mut ss = 0.0_f64;
113                for w in random.windows(2) {
114                    let a = y[c * n_rows + units[w[0]]] / expected_hits[w[0]];
115                    let b = y[c * n_rows + units[w[1]]] / expected_hits[w[1]];
116                    ss += (b - a) * (b - a);
117                }
118                draws as f64 / (2.0 * (draws as f64 - 1.0)) * ss
119            }
120        })
121        .collect();
122
123    Ok(Estimate {
124        total,
125        random_draws: draws,
126        sd_variance,
127    })
128}
chromy_core/estimate.rs

chromy_core/
estimate.rs
