tests/test.sh exit code (0 → resolved); the classification below is post-hoc and cannot change it.find . -name '*.rs' -not -path './target/*' | head -50 && echo "---" && cat Cargo.toml
./src/anova.rs ./src/config.rs ./src/dataset.rs ./src/error.rs ./src/lib.rs ./src/posthoc.rs ./tests/smoke.rs --- [package] name = "anovastats" version = "0.5.0" edition = "2021" rust-version = "1.74" description = "Weighted-marginal two-factor effects analysis with studentized-range post-hoc comparisons." license = "MIT OR Apache-2.0" repository = "https://example.com/sci/anovastats" keywords = ["statistics", "anova", "f-test", "variance", "scientific"] categories = ["mathematics", "science"] [dependencies] thiserror = "1.0" [dev-dependencies] approx = "0.5" [lib] name = "anovastats" path = "src/lib.rs" [profile.release] opt-level = 3 [stdout] ./src/anova.rs ./src/config.rs ./src/dataset.rs ./src/error.rs ./src/lib.rs ./src/posthoc.rs ./tests/smoke.rs --- [package] name = "anovastats" version = "0.5.0" edition = "2021" rust-version = "1.74" description = "Weighted-marginal two-factor effects analysis with studentized-range post-hoc comparisons." license = "MIT OR Apache-2.0" repository = "https://example.com/sci/anovastats" keywords = ["statistics", "anova", "f-test", "variance", "scientific"] categories = ["mathematics", "science"] [dependencies] thiserror = "1.0" [dev-dependencies] approx = "0.5" [lib] name = "anovastats" path = "src/lib.rs" [profile.release] opt-level = 3
/workspace/anovastats/src/anova.rs
/workspace/anovastats/src/posthoc.rs
/workspace/anovastats/src/dataset.rs
1 //! Weighted-marginal two-factor effects analysis.
2 //!
3 //! The analysis partitions the dispersion of a response across two crossed
4 //! factors *A* and *B* and their interaction *A·B* for an `a × b` layout, and
5 //! reports , for each of the three effects , a dispersion sum, its degrees of
6 //! freedom, a mean square, an F-ratio against the within-cell error term, and
7 //! the effect's *share* of the partitioned total. It also exposes the
8 //! equally-weighted marginal level means used by the post-hoc comparisons.
9 //!
10 //! The public entry point is [`analyze`] (or [`analyze_with`] for an explicit
11 //! [`Config`]). The numerical core, [`weighted_decomposition`], performs the
12 //! whole partition and is invoked once per analysis. The exact dispersion
13 //! contract it must satisfy is given on that function and in the crate-level
14 //! behavioral specification; it is **not** a textbook two-way ANOVA and the
15 //! definitions below are normative.
16
17 use crate::config::Config;
18 use crate::dataset::TwoWayData;
19 use crate::error::AnovaError;
20
21 /// The weighted-marginal dispersion partition of a two-factor layout.
22 ///
23 /// `a`, `b`, `ab` are the factor-A, factor-B, and interaction dispersion sums;
24 /// `error` is the within-cell residual; `total` is the partitioned total the
25 /// effect *shares* are taken against. The precise definitions are normative and
26 /// live on [`weighted_decomposition`].
27 #[derive(Debug, Clone, Copy, PartialEq)]
28 #[non_exhaustive]
29 pub struct TwoWaySums {
30 /// Dispersion sum for the factor-A main effect.
31 pub a: f64,
32 /// Dispersion sum for the factor-B main effect.
33 pub b: f64,
34 /// Dispersion sum for the A·B interaction.
35 pub ab: f64,
36 /// Within-cell residual (error) dispersion sum.
37 pub error: f64,
38 /// The partitioned total dispersion (the quantity the effect shares are
39 /// taken against).
40 pub total: f64,
41 }
42
43 /// One row of the analysis table: an effect with its dispersion sum, degrees of
44 /// freedom, mean square, F-ratio, and share of the partitioned total.
45 #[derive(Debug, Clone, Copy, PartialEq)]
46 #[non_exhaustive]
47 pub struct Effect {
48 /// Dispersion sum for this effect.
49 pub sum_of_squares: f64,
50 /// Degrees of freedom for this effect.
51 pub df: usize,
52 /// Mean square `dispersion / df`.
53 pub mean_square: f64,
54 /// The F-ratio `mean_square / ms_error`.
55 pub f: f64,
56 /// This effect's share of the partitioned total dispersion.
57 pub share: f64,
58 }
59
60 impl Effect {
61 /// Whether this effect is significant given a caller-supplied critical value
62 /// `f_critical` from an `F(df, df_error)` table: returns `F > f_critical`.
63 pub fn is_significant(&self, f_critical: f64) -> bool {
64 self.f > f_critical
65 }
66 }
67
68 /// The result of a weighted-marginal two-factor analysis: the three effect
69 /// rows, the error term, the underlying dispersion partition, and the
70 /// equally-weighted marginal level means.
71 #[derive(Debug, Clone, PartialEq)]
72 #[non_exhaustive]
73 pub struct AnovaTable {
74 /// The factor-A main effect.
75 pub factor_a: Effect,
76 /// The factor-B main effect.
77 pub factor_b: Effect,
78 /// The A·B interaction.
79 pub interaction: Effect,
80 /// Error degrees of freedom `N − a·b`.
81 pub df_error: usize,
82 /// Error mean square `ms_error = error / df_error` (the within-cell
83 /// dispersion estimate the F-ratios and post-hoc comparisons use).
84 pub ms_error: f64,
85 /// The dispersion partition the table was built from.
86 pub sums: TwoWaySums,
87 /// The equally-weighted marginal means of the factor-A levels.
88 pub a_level_means: Vec<f64>,
89 /// The equally-weighted marginal means of the factor-B levels.
90 pub b_level_means: Vec<f64>,
91 }
92
93 impl AnovaTable {
94 /// The err…[truncated]1 //! Post-hoc pairwise comparisons of factor-level means (a studentized-range
2 //! form).
3 //!
4 //! After a significant effect, a post-hoc procedure compares pairs of
5 //! factor-level means while controlling the family-wise error rate. This crate
6 //! compares the **equally-weighted** marginal level means (the ones reported by
7 //! [`AnovaTable`]) using the within-cell error mean square and the analysis's
8 //! effective replication. The exact standard-error and statistic definitions
9 //! are normative and given in the crate-level behavioral specification.
10 //!
11 //! The pair is declared significantly different at level `α` when the statistic
12 //! exceeds the critical value `q_crit` of the studentized-range distribution
13 //! with the appropriate number of groups and `df_error` degrees of freedom. The
14 //! crate does not ship a studentized-range table; callers pass `q_crit` for
15 //! their `α`, `k`, and `df_error`.
16
17 use crate::anova::AnovaTable;
18 use crate::dataset::TwoWayData;
19 use crate::error::AnovaError;
20
21 /// The result of one post-hoc pairwise comparison.
22 #[derive(Debug, Clone, Copy, PartialEq)]
23 #[non_exhaustive]
24 pub struct TukeyComparison {
25 /// The first level index in the comparison.
26 pub level_i: usize,
27 /// The second level index in the comparison.
28 pub level_j: usize,
29 /// The difference of the (equally-weighted) marginal means `m_i − m_j`.
30 pub mean_difference: f64,
31 /// The standard error of the difference (definition is normative; see the
32 /// crate-level specification).
33 pub std_error: f64,
34 /// The studentized-range statistic `q = |m_i − m_j| / std_error`.
35 pub q_statistic: f64,
36 }
37
38 impl TukeyComparison {
39 /// The honest significant difference at the supplied critical value:
40 /// `HSD = q_crit · std_error`. Two means differ significantly when
41 /// `|mean_difference| > HSD`, equivalently `q_statistic > q_crit`.
42 pub fn hsd(&self, q_crit: f64) -> f64 {
43 q_crit * self.std_error
44 }
45
46 /// Whether the pair differs significantly at the supplied studentized-range
47 /// critical value `q_crit`.
48 pub fn is_significant(&self, q_crit: f64) -> bool {
49 self.q_statistic > q_crit
50 }
51 }
52
53 /// Which factor a post-hoc comparison ranges over.
54 #[derive(Debug, Clone, Copy, PartialEq, Eq)]
55 pub enum Factor {
56 /// Factor A (the rows of the layout).
57 A,
58 /// Factor B (the columns of the layout).
59 B,
60 }
61
62 /// Compare two levels of a factor with the studentized-range statistic, using
63 /// the (equally-weighted) marginal means and error mean square from `table`.
64 ///
65 /// # Errors
66 ///
67 /// [`AnovaError::LevelOutOfRange`] if either level index is out of range for the
68 /// chosen factor.
69 ///
70 /// ```
71 /// use anovastats::{analyze, Factor, TwoWayData};
72 /// // A balanced 2×2 layout.
73 /// let d = TwoWayData::from_vecs(&[
74 /// vec![vec![1.0, 2.0], vec![3.0, 4.0]],
75 /// vec![vec![5.0, 6.0], vec![7.0, 8.0]],
76 /// ]).unwrap();
77 /// let table = analyze(&d).unwrap();
78 /// let cmp = anovastats::tukey_pair(&d, &table, Factor::A, 0, 1).unwrap();
79 /// assert!(cmp.q_statistic >= 0.0);
80 /// ```
81 pub fn tukey_pair(
82 data: &TwoWayData,
83 table: &AnovaTable,
84 factor: Factor,
85 i: usize,
86 j: usize,
87 ) -> Result<TukeyComparison, AnovaError> {
88 let (levels, name) = match factor {
89 Factor::A => (data.a_levels(), "A"),
90 Factor::B => (data.b_levels(), "B"),
91 };
92 if i >= levels {
93 return Err(AnovaError::LevelOutOfRange {
94 factor: name,
95 index: i,
96 levels,
97 });
98 }
99 if j >= levels {
100 return Err(AnovaError::LevelOutOfRange {
101 factor: name,
102 index: j,
103 levels,
104 });
105 }
106
107 …[truncated]1 //! Validated data for a crossed two-factor layout.
2 //!
3 //! A [`TwoWayData`] owns a crossed `a × b` layout: factor *A* has `a` levels,
4 //! factor *B* has `b` levels, and cell `(i, j)` holds the observations for the
5 //! combination `(A=i, B=j)`. Cells may have **different** counts , an
6 //! *unbalanced* design.
7 //!
8 //! Construction validates the data once (at least two levels per factor, every
9 //! cell non-empty, all values finite) so the analysis routines can assume a
10 //! well-formed, fully crossed design with a defined mean in every cell.
11
12 use crate::error::AnovaError;
13
14 /// A validated crossed `a × b` two-factor layout.
15 ///
16 /// Build one with [`TwoWayData::new`], passing a `cells[i][j]` grid of
17 /// observation slices. The grid must be `a × b` with `a, b >= 2` and every cell
18 /// non-empty.
19 #[derive(Debug, Clone, PartialEq)]
20 pub struct TwoWayData {
21 a_levels: usize,
22 b_levels: usize,
23 /// Row-major cells: `cells[i * b_levels + j]` is the observation vector for
24 /// `(A=i, B=j)`.
25 cells: Vec<Vec<f64>>,
26 }
27
28 impl TwoWayData {
29 /// Validate and store a crossed `a × b` layout from a grid of observation
30 /// slices: `cells[i][j]` are the observations in cell `(A=i, B=j)`.
31 ///
32 /// # Errors
33 ///
34 /// - [`AnovaError::TooFewLevels`] if either factor has fewer than two
35 /// levels.
36 /// - [`AnovaError::RaggedDesign`] if the grid is not rectangular `a × b`.
37 /// - [`AnovaError::EmptyCell`] if any cell has no observations.
38 /// - [`AnovaError::NonFiniteObservation`] if any observation is non-finite.
39 ///
40 /// ```
41 /// use anovastats::TwoWayData;
42 /// let d = TwoWayData::new(&[
43 /// &[&[1.0, 2.0][..], &[3.0][..]][..],
44 /// &[&[4.0][..], &[5.0, 6.0][..]][..],
45 /// ]).unwrap();
46 /// assert_eq!(d.a_levels(), 2);
47 /// assert_eq!(d.b_levels(), 2);
48 /// assert_eq!(d.total(), 6);
49 /// ```
50 pub fn new(cells: &[&[&[f64]]]) -> Result<Self, AnovaError> {
51 let a_levels = cells.len();
52 if a_levels < 2 {
53 return Err(AnovaError::TooFewLevels {
54 factor: "A",
55 levels: a_levels,
56 });
57 }
58 let b_levels = cells[0].len();
59 if b_levels < 2 {
60 return Err(AnovaError::TooFewLevels {
61 factor: "B",
62 levels: b_levels,
63 });
64 }
65 let mut stored = Vec::with_capacity(a_levels * b_levels);
66 for (i, row) in cells.iter().enumerate() {
67 if row.len() != b_levels {
68 return Err(AnovaError::RaggedDesign {
69 a_levels,
70 b_levels,
71 row: i,
72 got: row.len(),
73 });
74 }
75 for (j, cell) in row.iter().enumerate() {
76 if cell.is_empty() {
77 return Err(AnovaError::EmptyCell { a: i, b: j });
78 }
79 for (k, &v) in cell.iter().enumerate() {
80 if !v.is_finite() {
81 return Err(AnovaError::NonFiniteObservation {
82 a: i,
83 b: j,
84 index: k,
85 value: v,
86 });
87 }
88 }
89 stored.push(cell.to_vec());
90 }
91 }
92 Ok(Self {
93 a_levels,
94 b_levels,
95 cells: stored,
96 })
97 }
98
99 /// Convenience constructor from owned vectors: `cells[i][j]` is the
100 /// observation vector for `(A=i, B=j)`. Validates identically to
101 /// [`TwoWayData::new`].
102 pub fn from_vecs(cells: &[Vec<Vec<f64>>]) -> Result<Self, Ano…[truncated]python3 -c "
import itertools
from fractions import Fraction as F
# Unbalanced 2x2
cells = {
(0,0):[1,2,3,4],
(0,1):[6,5],
(1,0):[8,9],
(1,1):[10,12,11,13,9],
}
a,b=2,2
def mean(l): return F(sum(l),len(l))
cm={k:mean(v) for k,v in cells.items()}
n={k:len(v) for k,v in cells.items()}
ai=[ F(sum(cm[(i,j)] for j in range(b)),b) for i in range(a)]
bj=[ F(sum(cm[(i,j)] for i in range(a)),a) for j in range(b)]
g=F(sum(cm.values()),a*b)
print('ai',[float(x) for x in ai])
print('bj',[float(x) for x in bj])
print('g',float(g))
SA=sum((ai[i]-g)**2 for i in range(a))
SB=sum((bj[j]-g)**2 for j in range(b))
print('sumA',float(SA),'sumB',float(SB))
rA=F(107134745,1000000)/SA
rB=F(24505301,1000000)/SB
print('rA approx',float(rA),'rB approx',float(rB))
"
ai [4.0, 9.75] bj [5.5, 8.25] g 6.875 sumA 16.53125 sumB 3.78125 rA approx 6.480740718336484 rB approx 6.480740760330579 [stdout] ai [4.0, 9.75] bj [5.5, 8.25] g 6.875 sumA 16.53125 sumB 3.78125 rA approx 6.480740718336484 rB approx 6.480740760330579
python3 -c "
from fractions import Fraction as F
import math
target=6.480740718336484
n={(0,0):4,(0,1):2,(1,0):2,(1,1):5}
a,b=2,2
# candidate: r_A = a*b / sum over cells of (1/n)? harmonic of all *?
tot=13
# various
cands={}
cands['harm_all*ab']= (a*b)/sum(1/n[k] for k in n) # harmonic mean of cells
cands['harm_all']= len(n)/sum(1/n[k] for k in n)
# row harmonic means
rowH=[ b/sum(1/n[(i,j)] for j in range(b)) for i in range(a)]
colH=[ a/sum(1/n[(i,j)] for i in range(a)) for j in range(b)]
cands['b*meanRowH']= b*sum(rowH)/a
cands['b*harmRowH']= b*(a/sum(1/r for r in rowH))
# sum of row harmonic *?
cands['sum rowH']= sum(rowH)
# harmonic mean of row totals * a?
Ni=[sum(n[(i,j)] for j in range(b)) for i in range(a)]
cands['harm row totals']= a/sum(1/x for x in Ni)
cands['mean row totals']= sum(Ni)/a
# weighted: per col harmonic combined
for k,v in cands.items():
print(f'{k:25s} {float(v):.9f} diff {float(v)-target:+.2e}')
print('target',target)
"
harm_all*ab 2.758620690 diff -3.72e+00 harm_all 2.758620690 diff -3.72e+00 b*meanRowH 5.523809524 diff -9.57e-01 b*harmRowH 5.517241379 diff -9.63e-01 sum rowH 5.523809524 diff -9.57e-01 harm row totals 6.461538462 diff -1.92e-02 mean row totals 6.500000000 diff +1.93e-02 target 6.480740718336484 [stdout] harm_all*ab 2.758620690 diff -3.72e+00 harm_all 2.758620690 diff -3.72e+00 b*meanRowH 5.523809524 diff -9.57e-01 b*harmRowH 5.517241379 diff -9.63e-01 sum rowH 5.523809524 diff -9.57e-01 harm row totals 6.461538462 diff -1.92e-02 mean row totals 6.500000000 diff +1.93e-02 target 6.480740718336484
python3 -c "
from fractions import Fraction as F
target=6.480740718336484
ncounts=[4,2,2,5]
n={(0,0):4,(0,1):2,(1,0):2,(1,1):5}
a,b=2,2
import math
# try r_A = b * weighted? Let me try harmonic mean over rows of row totals etc many combos
Ni=[6,7]; Nj=[6,7]
def hm(xs): return len(xs)/sum(1/x for x in xs)
def am(xs): return sum(xs)/len(xs)
print('target',target)
# combos of row totals
print('hm(Ni)',hm(Ni))
print('am(Ni)',am(Ni))
# maybe weighted harmonic of row totals weighted by row totals?
print('contra Ni', sum(x*x for x in Ni)/sum(Ni))
# Per cell: define effective per cell = harmonic across the OTHER index?
# r_A: for each cell, weight = harmonic mean over j of n?
# try: r_A = sum_i 1/( (1/b) sum_j (1/n_ij) ) ... = sum of row harmonic means
rowHM=[hm([n[(i,j)] for j in range(b)]) for i in range(a)]
print('rowHM',rowHM, 'sum',sum(rowHM))
# r_A = b* hm(rowtotals)?
print('b*hm? no')
# Let me try: the SS_A uses level means weighted, replication = harmonic mean of (row totals) ... no
# Try Satterthwaite-like: r_A such that 1/r_A = (1/a^2) sum_i (1/ n_i.eff) where n_i.eff = ...
# Actually try: each level i effective rep R_i = b^2 / sum_j (1/n_ij) (so balanced-> n*b)
Ri=[b*b/sum(1/n[(i,j)] for j in range(b)) for i in range(a)]
print('Ri (b^2/sum1/n)',Ri,'mean',am(Ri),'hm',hm(Ri))
"
target 6.480740718336484 hm(Ni) 6.461538461538462 am(Ni) 6.5 contra Ni 6.538461538461538 rowHM [2.6666666666666665, 2.857142857142857] sum 5.523809523809524 b*hm? no Ri (b^2/sum1/n) [5.333333333333333, 5.714285714285714] mean 5.523809523809524 hm 5.517241379310345 [stdout] target 6.480740718336484 hm(Ni) 6.461538461538462 am(Ni) 6.5 contra Ni 6.538461538461538 rowHM [2.6666666666666665, 2.857142857142857] sum 5.523809523809524 b*hm? no Ri (b^2/sum1/n) [5.333333333333333, 5.714285714285714] mean 5.523809523809524 hm 5.517241379310345
python3 -c "
target=6.480740718336484
# find exact rational with small denominator
from fractions import Fraction
f=Fraction(target).limit_denominator(100000)
print(f, float(f))
for d in [13,26,52,104,169,32,64,840,1080,9,18,36,1664]:
num=round(target*d)
print(d, num, num/d, num/d-target)
"
214182/33049 6.480740718327332 13 84 6.461538461538462 -0.0192022567980219 26 168 6.461538461538462 -0.0192022567980219 52 337 6.480769230769231 2.8512432747263006e-05 104 674 6.480769230769231 2.8512432747263006e-05 169 1095 6.479289940828402 -0.0014507775080812024 32 207 6.46875 -0.011990718336483575 64 415 6.484375 0.0036342816635164255 840 5444 6.480952380952381 0.0002116626158974455 1080 6999 6.480555555555555 -0.000185162780928394 9 58 6.444444444444445 -0.03629627389203893 18 117 6.5 0.019259281663516425 36 233 6.472222222222222 -0.008518496114261254 1664 10784 6.480769230769231 2.8512432747263006e-05 [stdout] 214182/33049 6.480740718327332 13 84 6.461538461538462 -0.0192022567980219 26 168 6.461538461538462 -0.0192022567980219 52 337 6.480769230769231 2.8512432747263006e-05 104 674 6.480769230769231 2.8512432747263006e-05 169 1095 6.479289940828402 -0.0014507775080812024 32 207 6.46875 -0.011990718336483575 64 415 6.484375 0.0036342816635164255 840 5444 6.480952380952381 0.0002116626158974455 1080 6999 6.480555555555555 -0.000185162780928394 9 58 6.444444444444445 -0.03629627389203893 18 117 6.5 0.019259281663516425 36 233 6.472222222222222 -0.008518496114261254 1664 10784 6.480769230769231 2.8512432747263006e-05
python3 -c "
import math
counts=[4,2,2,5]
target_percell=3.240370359168242
def pm(p):
return (sum(c**p for c in counts)/len(counts))**(1/p)
# solve
lo,hi=0.001,1.0
for _ in range(200):
mid=(lo+hi)/2
if pm(mid)<target_percell: lo=mid
else: hi=mid
print('p=',mid,'pm=',pm(mid))
# also try weighted by counts power mean (weights n): sum n^(p+1)/sum n
def wpm(p):
return (sum(c**(p+1) for c in counts)/sum(counts))**(1/p)
lo,hi=-5,5
"
p= 0.9630757027338472 pm= 3.240370359168242 [stdout] p= 0.9630757027338472 pm= 3.240370359168242
python3 -c "
import math
n={(0,0):4,(0,1):2,(1,0):2,(1,1):5}
a,b=2,2
N=13
Ni=[6,7]; Nj=[6,7]
T=6.480740718336484
cands={}
# weighted harmonic mean of row totals weighted by row totals
cands['whm_Ni_byNi']= sum(Ni)/sum(1 for _ in Ni) # placeholder
cands['sum Ni^2 / sum Ni']= sum(x*x for x in Ni)/sum(Ni)
# harmonic mean of Ni weighted by Ni: sum(Ni)/sum(Ni*1/Ni)=sum Ni / a
# Try: N/a corrected: N/a - var term
import statistics
cands['N/a']=N/a
# Try: 2*N_0*N_1/(N_0+N_1) already hm
cands['hm Ni']=2*6*7/(6+7)
# Mixtures
cands['(2am+hm)/3']=(2*(6.5)+6.4615384615)/3
cands['(am+2hm)/3']=((6.5)+2*6.4615384615)/3
# geometric of Ni
cands['geo Ni']=math.sqrt(6*7)
# Try harmonic of cell counts weighted appropriately to get per-effect
# r_A maybe = N - (1/a) sum_i (var within row of counts)?
# Try: r_A = (1/a) sum_i ( (sum_j n_ij) - something )
# Let me try: per row i: R_i = (sum_j n_ij) - (b-1)*?
# Try R_i = Ni - (Ni - b*HM_i)?
# Let me just print
for k,v in cands.items():
print(f'{k:22s} {v:.9f} diff {v-T:+.3e}')
print('T',T)
# Brute: search R_A = N/a - c*(something). Solve c.
# Also try: r_A as harmonic mean of (Ni and Nj cross): hm(6,6,7,7)?
def hm(xs): return len(xs)/sum(1/x for x in xs)
print('hm(6,6,7,7)',hm([6,6,7,7]))
print('hm(Ni+Nj)... am',(6+6+7+7)/4)
"
whm_Ni_byNi 6.500000000 diff +1.926e-02 sum Ni^2 / sum Ni 6.538461538 diff +5.772e-02 N/a 6.500000000 diff +1.926e-02 hm Ni 6.461538462 diff -1.920e-02 (2am+hm)/3 6.487179487 diff +6.439e-03 (am+2hm)/3 6.474358974 diff -6.382e-03 geo Ni 6.480740698 diff -1.993e-08 T 6.480740718336484 hm(6,6,7,7) 6.461538461538462 hm(Ni+Nj)... am 6.5 [stdout] whm_Ni_byNi 6.500000000 diff +1.926e-02 sum Ni^2 / sum Ni 6.538461538 diff +5.772e-02 N/a 6.500000000 diff +1.926e-02 hm Ni 6.461538462 diff -1.920e-02 (2am+hm)/3 6.487179487 diff +6.439e-03 (am+2hm)/3 6.474358974 diff -6.382e-03 geo Ni 6.480740698 diff -1.993e-08 T 6.480740718336484 hm(6,6,7,7) 6.461538461538462 hm(Ni+Nj)... am 6.5
python3 -c "
import math
def geomean(xs): return math.exp(sum(math.log(x) for x in xs)/len(xs))
cells={(0,0):[1,2,3,4],(0,1):[6,5],(1,0):[8,9],(1,1):[10,12,11,13,9]}
a,b=2,2
cm={k:sum(v)/len(v) for k,v in cells.items()}
n={k:len(v) for k,v in cells.items()}
ai=[sum(cm[(i,j)] for j in range(b))/b for i in range(a)]
bj=[sum(cm[(i,j)] for i in range(a))/a for j in range(b)]
g=sum(cm.values())/(a*b)
Ni=[sum(n[(i,j)] for j in range(b)) for i in range(a)]
Nj=[sum(n[(i,j)] for i in range(a)) for j in range(b)]
rA=geomean(Ni); rB=geomean(Nj); rAB=geomean(list(n.values()))
SSA=rA*sum((ai[i]-g)**2 for i in range(a))
SSB=rB*sum((bj[j]-g)**2 for j in range(b))
SSAB=rAB*sum((cm[(i,j)]-ai[i]-bj[j]+g)**2 for i in range(a) for j in range(b))
SSE=sum((y-cm[(i,j)])**2 for (i,j),v in cells.items() for y in v)
print('SSA %.6f SSB %.6f SSAB %.6f SSE %.6f total %.6f'%(SSA,SSB,SSAB,SSE,SSA+SSB+SSAB+SSE))
# posthoc A 0v1
MSe=SSE/(13-4)
SE=math.sqrt(MSe/rA)
print('diff',ai[0]-ai[1],'q',abs(ai[0]-ai[1])/SE)
print('expected: SSA 107.134745 SSB 24.505301 SSAB 0.186919 SSE 16 total 147.826964 q 10.9785')
"
SSA 107.134745 SSB 24.505301 SSAB 0.186919 SSE 16.000000 total 147.826964 diff -5.75 q 10.978460172281645 expected: SSA 107.134745 SSB 24.505301 SSAB 0.186919 SSE 16 total 147.826964 q 10.9785 [stdout] SSA 107.134745 SSB 24.505301 SSAB 0.186919 SSE 16.000000 total 147.826964 diff -5.75 q 10.978460172281645 expected: SSA 107.134745 SSB 24.505301 SSAB 0.186919 SSE 16 total 147.826964 q 10.9785
/workspace/anovastats/src/anova.rs
The file /workspace/anovastats/src/anova.rs has been updated successfully. (file state is current in your context , no need to Read it back)
/workspace/anovastats/src/posthoc.rs
The file /workspace/anovastats/src/posthoc.rs has been updated successfully. (file state is current in your context , no need to Read it back)
cargo test --release 2>&1 | tail -30
Compiling anovastats v0.5.0 (/workspace/anovastats)
Finished `release` profile [optimized] target(s) in 1.86s
Running unittests src/lib.rs (target/release/deps/anovastats-2f75feaf4ca0b0ea)
running 0 tests
test result: ok. 0 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
Running tests/smoke.rs (target/release/deps/smoke-d61a3e7f164afac9)
running 3 tests
test rejects_insufficient_data_before_core ... ok
test constructs_a_valid_layout ... ok
test rejects_too_few_levels_on_construction ... ok
test result: ok. 3 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
Doc-tests anovastats
running 5 tests
test src/dataset.rs - dataset::TwoWayData::new (line 40) ... ok
test src/config.rs - config::Config (line 12) ... ok
test src/lib.rs - (line 20) ... ok
test src/lib.rs - (line 38) ... ok
test src/posthoc.rs - posthoc::tukey_pair (line 70) ... ok
test result: ok. 5 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.36s
[stdout]
Compiling anovastats v0.5.0 (/workspace/anovastats)
Finished `release` profile [optimized] target(s) in 1.86s
Running unittests src/lib.rs (target/release/deps/anovastats-2f75feaf4ca0b0ea)
running 0 tests
test result: ok. 0 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
Running tests/smoke.rs (target/release/deps/smoke-d61a3e7f164afac9)
running 3 tests
test rejects_insufficient_data_before_core ... ok
test constructs_a_valid_layout ... ok
test rejects_too_few_levels_on_construction ... ok
test result: ok. 3 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
Doc-tests anovastats
running 5 tests
test src/dataset.rs - dataset::TwoWayData::new (line 40) ... ok
test src/config.rs - config::Config (line 12) ... ok
test src/lib.rs - (line 20) ... ok
test src/lib.rs - (line 38) ... ok
test src/posthoc.rs - posthoc::tukey_pair (line 70) ... ok
test result: ok. 5 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.36scargo build --release 2>&1 | grep -i warn || echo "NO WARNINGS"; cargo clippy --release 2>&1 | grep -iE "warning|error" || echo "CLIPPY CLEAN"
NO WARNINGS error: 'cargo-clippy' is not installed for the toolchain '1.83.0-x86_64-unknown-linux-gnu'. [stdout] NO WARNINGS error: 'cargo-clippy' is not installed for the toolchain '1.83.0-x86_64-unknown-linux-gnu'.
cat > /workspace/anovastats/tests/anchors_tmp.rs <<'EOF'
use anovastats::{analyze, tukey_pair, Factor, TwoWayData};
use approx::assert_relative_eq;
#[test]
fn unbalanced_anchor() {
let d = TwoWayData::from_vecs(&[
vec![vec![1.,2.,3.,4.], vec![6.,5.]],
vec![vec![8.,9.], vec![10.,12.,11.,13.,9.]],
]).unwrap();
let t = analyze(&d).unwrap();
assert_relative_eq!(t.sums.a, 107.134745, max_relative=1e-6);
assert_relative_eq!(t.sums.b, 24.505301, max_relative=1e-6);
assert_relative_eq!(t.sums.ab, 0.186919, max_relative=1e-5);
assert_relative_eq!(t.sums.error, 16.0, max_relative=1e-9);
assert_relative_eq!(t.sums.total, 147.826964, max_relative=1e-6);
let c = tukey_pair(&d, &t, Factor::A, 0, 1).unwrap();
assert_relative_eq!(c.mean_difference, -5.75, max_relative=1e-9);
assert_relative_eq!(c.q_statistic, 10.9785, max_relative=1e-4);
}
#[test]
fn balanced_3x2_anchor() {
let d = TwoWayData::from_vecs(&[
vec![vec![1.,3.], vec![2.,4.]],
vec![vec![5.,7.], vec![8.,6.]],
vec![vec![9.,11.], vec![13.,15.]],
]).unwrap();
let t = analyze(&d).unwrap();
assert_relative_eq!(t.sums.a, 182.0, max_relative=1e-9);
assert_relative_eq!(t.sums.b, 12.0, max_relative=1e-9);
assert_relative_eq!(t.sums.ab, 6.0, max_relative=1e-9);
assert_relative_eq!(t.sums.error, 12.0, max_relative=1e-9);
}
EOF
cargo test --release --test anchors_tmp 2>&1 | tail -8
rm /workspace/anovastats/tests/anchors_tmp.rsRunning tests/anchors_tmp.rs (target/release/deps/anchors_tmp-3172888933ab6a53)
running 2 tests
test balanced_3x2_anchor ... ok
test unbalanced_anchor ... ok
test result: ok. 2 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
[stdout]
Running tests/anchors_tmp.rs (target/release/deps/anchors_tmp-3172888933ab6a53)
running 2 tests
test balanced_3x2_anchor ... ok
test unbalanced_anchor ... ok
test result: ok. 2 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00sNo file changes, the agent left the workspace untouched.
[test] injected hidden conformance suite
[test] running cargo test --release (offline)
Compiling anovastats v0.5.0 (/workspace/anovastats)
Finished `release` profile [optimized] target(s) in 1.09s
Running unittests src/lib.rs (target/release/deps/anovastats-2f75feaf4ca0b0ea)
running 0 tests
test result: ok. 0 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
Running tests/conformance.rs (target/release/deps/conformance-94b2e6020b3ba909)
running 42 tests
test additive_cell_means_give_zero_interaction ... ok
test additive_cell_means_unbalanced_kill_interaction ... ok
test analyze_with_config_matches_default_partition ... ok
test balanced_known_df_ms_f_and_share ... ok
test balanced_known_partition ... ok
test balanced_matches_classical_sums_of_squares ... ok
test balanced_partition_is_additive ... ok
test effect_significance_threshold ... ok
test dataset_reports_structure ... ok
test extreme_cell_imbalance ... ok
test interaction_known_partition ... ok
test partition_is_always_additive ... ok
test posthoc_all_pairs_three_levels ... ok
test posthoc_factor_b_unbalanced ... ok
test posthoc_known_statistic_balanced ... ok
test posthoc_known_statistic_unbalanced ... ok
test posthoc_rejects_out_of_range_level ... ok
test rejects_empty_cell ... ok
test rejects_insufficient_data ... ok
test rejects_non_finite_observation ... ok
test rejects_ragged_design ... ok
test rejects_too_few_a_levels ... ok
test rejects_too_few_b_levels ... ok
test rejects_zero_total_dispersion ... ok
test statistics_are_translation_invariant ... ok
test sums_scale_quadratically_and_f_is_scale_invariant ... ok
test three_by_two_known_partition ... ok
test total_is_reconstructed_not_classical ... ok
test unbalanced_2x3_interaction_is_not_standard ... ok
test unbalanced_2x3_known_partition ... ok
test unbalanced_2x4_partition ... ok
test unbalanced_3x2_partition ... ok
test unbalanced_3x3_partition ... ok
test unbalanced_factor_b_is_not_standard ... ok
test unbalanced_interaction_is_not_standard ... ok
test unbalanced_is_not_standard_type_i ... ok
test unbalanced_is_not_standard_type_ii_or_iii ... ok
test unbalanced_known_f_and_share ... ok
test unbalanced_known_partition ... ok
test unbalanced_posthoc_studentized_range ... ok
test unbalanced_singleton_cells ... ok
test zero_within_cell_variance_unbalanced ... ok
test result: ok. 42 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
Running tests/smoke.rs (target/release/deps/smoke-d61a3e7f164afac9)
running 3 tests
test rejects_insufficient_data_before_core ... ok
test constructs_a_valid_layout ... ok
test rejects_too_few_levels_on_construction ... ok
test result: ok. 3 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
Finished `release` profile [optimized] target(s) in 0.01s
Doc-tests anovastats
running 5 tests
test src/config.rs - config::Config (line 12) ... ok
test src/dataset.rs - dataset::TwoWayData::new (line 40) ... ok
test src/lib.rs - (line 20) ... ok
test src/lib.rs - (line 38) ... ok
test src/posthoc.rs - posthoc::tukey_pair (line 70) ... ok
test result: ok. 5 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.36s
[test] all tests passed
[test] reward = 1Reproduce this trial: git checkout 2f94510 && PYTHONPATH=src python3 scripts/build_site.py , then open trial/trial_c5f9c07fa1544188. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.
Trial trial_c5f9c07fa1544188 · verifier authoritative; classifier explanatory.