diff --git a/.gitignore b/.gitignore
index 573ede2..40a3f94 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,5 @@
 /target
 /Cargo.lock
 
+.justfile
 *.svg
diff --git a/Cargo.toml b/Cargo.toml
index d579204..496042c 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -3,10 +3,18 @@ name = "trueskill-tt"
 version = "0.1.0"
 edition = "2021"
 
+[lib]
+bench = false
+
+[[bench]]
+name = "batch"
+harness = false
+
 [dependencies]
 approx = { version = "0.5.1", optional = true }
 
 [dev-dependencies]
+criterion = "0.4"
 plotters = { version = "0.3.4", default-features = false, features = ["svg_backend", "all_elements", "all_series"] }
 plotters-backend = "0.3.4"
 time = { version = "0.3.17", features = ["parsing"] }
diff --git a/README.md b/README.md
index 7672b0f..dddd7f8 100644
--- a/README.md
+++ b/README.md
@@ -13,6 +13,7 @@ Rust port of [TrueSkillThroughTime.py](https://github.com/glandfried/TrueSkillTh
 
 - [x] Implement approx for Gaussian
 - [x] Add more tests from `TrueSkillThroughTime.jl`
+- [ ] Benchmark Batch::iteration()
 - [ ] Time needs to be an enum so we can have multiple states (see `batch::compute_elapsed()`)
 - [ ] Add examples (use same TrueSkillThroughTime.(py|jl))
 - [ ] Add Observer (see [argmin](https://docs.rs/argmin/latest/argmin/core/trait.Observe.html) for inspiration)
diff --git a/benches/batch.rs b/benches/batch.rs
new file mode 100644
index 0000000..5dc4924
--- /dev/null
+++ b/benches/batch.rs
@@ -0,0 +1,62 @@
+use std::collections::HashMap;
+
+use criterion::{criterion_group, criterion_main, Criterion};
+use trueskill_tt::{
+    agent::Agent, batch::Batch, gaussian::Gaussian, player::Player, IndexMap, BETA, GAMMA, MU,
+    P_DRAW, SIGMA,
+};
+
+fn criterion_benchmark(criterion: &mut Criterion) {
+    let mut index = IndexMap::new();
+
+    let a = index.get_or_create("a");
+    let b = index.get_or_create("b");
+    let c = index.get_or_create("c");
+
+    let agents = {
+        let mut map = HashMap::new();
+
+        map.insert(
+            a,
+            Agent {
+                player: Player::new(Gaussian::new(MU, SIGMA), BETA, GAMMA),
+                ..Default::default()
+            },
+        );
+        map.insert(
+            b,
+            Agent {
+                player: Player::new(Gaussian::new(MU, SIGMA), BETA, GAMMA),
+                ..Default::default()
+            },
+        );
+        map.insert(
+            c,
+            Agent {
+                player: Player::new(Gaussian::new(MU, SIGMA), BETA, GAMMA),
+                ..Default::default()
+            },
+        );
+
+        map
+    };
+
+    let mut composition = Vec::new();
+    let mut results = Vec::new();
+    let mut weights = Vec::new();
+
+    for _ in 0..100 {
+        composition.push(vec![vec![a], vec![b]]);
+        results.push(vec![1.0, 0.0]);
+        weights.push(vec![vec![1.0], vec![1.0]]);
+    }
+
+    let mut batch = Batch::new(composition, results, weights, 1, P_DRAW, &agents);
+
+    criterion.bench_function("Batch::iteration", |b| {
+        b.iter(|| batch.iteration(0, &agents))
+    });
+}
+
+criterion_group!(benches, criterion_benchmark);
+criterion_main!(benches);
diff --git a/src/agent.rs b/src/agent.rs
index 671be40..c77cac6 100644
--- a/src/agent.rs
+++ b/src/agent.rs
@@ -1,10 +1,10 @@
 use crate::{gaussian::Gaussian, player::Player, N_INF};
 
 #[derive(Debug)]
-pub(crate) struct Agent {
-    pub(crate) player: Player,
-    pub(crate) message: Gaussian,
-    pub(crate) last_time: i64,
+pub struct Agent {
+    pub player: Player,
+    pub message: Gaussian,
+    pub last_time: i64,
 }
 
 impl Agent {
diff --git a/src/batch.rs b/src/batch.rs
index 9a0db6a..d974214 100644
--- a/src/batch.rs
+++ b/src/batch.rs
@@ -109,7 +109,7 @@ pub struct Batch {
 }
 
 impl Batch {
-    pub(crate) fn new(
+    pub fn new(
         composition: Vec<Vec<Vec<Index>>>,
         results: Vec<Vec<f64>>,
         weights: Vec<Vec<Vec<f64>>>,
@@ -269,13 +269,6 @@ impl Batch {
         self.iteration(from, agents);
     }
 
-    // TODO(anders): Use Item::posterior() instead.
-    pub fn posterior(&self, agent: Index) -> Gaussian {
-        let skill = &self.skills[&agent];
-
-        skill.likelihood * skill.backward * skill.forward
-    }
-
     pub(crate) fn posteriors(&self) -> HashMap<Index, Gaussian> {
         self.skills
             .iter()
@@ -283,7 +276,7 @@ impl Batch {
             .collect::<HashMap<_, _>>()
     }
 
-    pub(crate) fn iteration(&mut self, from: usize, agents: &HashMap<Index, Agent>) {
+    pub fn iteration(&mut self, from: usize, agents: &HashMap<Index, Agent>) {
         for event in self.events.iter_mut().skip(from) {
             let teams = event.within_priors(false, false, &self.skills, agents);
             let result = event.outputs();
diff --git a/src/history.rs b/src/history.rs
index 78da71b..9dddd70 100644
--- a/src/history.rs
+++ b/src/history.rs
@@ -460,7 +460,7 @@ mod tests {
 
         assert_ulps_eq!(observed, expected, epsilon = 0.000001);
 
-        let observed = h.batches[1].posterior(a);
+        let observed = h.batches[1].skills[&a].posterior();
 
         let p = Game::new(
             h.batches[1].events[0].within_priors(false, false, &h.batches[1].skills, &h.agents),
@@ -508,12 +508,12 @@ mod tests {
         h1.add_events_with_prior(composition, results, times, vec![], priors);
 
         assert_ulps_eq!(
-            h1.batches[0].posterior(a),
+            h1.batches[0].skills[&a].posterior(),
             Gaussian::new(22.904409, 6.010330),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h1.batches[0].posterior(c),
+            h1.batches[0].skills[&c].posterior(),
             Gaussian::new(25.110318, 5.866311),
             epsilon = 1e-6
         );
@@ -521,12 +521,12 @@ mod tests {
         h1.convergence(ITERATIONS, EPSILON, false);
 
         assert_ulps_eq!(
-            h1.batches[0].posterior(a),
+            h1.batches[0].skills[&a].posterior(),
             Gaussian::new(25.000000, 5.419212),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h1.batches[0].posterior(c),
+            h1.batches[0].skills[&c].posterior(),
             Gaussian::new(25.000000, 5.419212),
             epsilon = 1e-6
         );
@@ -553,12 +553,12 @@ mod tests {
         h2.add_events_with_prior(composition, results, times, vec![], priors);
 
         assert_ulps_eq!(
-            h2.batches[2].posterior(a),
+            h2.batches[2].skills[&a].posterior(),
             Gaussian::new(22.903522, 6.011017),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h2.batches[2].posterior(c),
+            h2.batches[2].skills[&c].posterior(),
             Gaussian::new(25.110702, 5.866811),
             epsilon = 1e-6
         );
@@ -566,12 +566,12 @@ mod tests {
         h2.convergence(ITERATIONS, EPSILON, false);
 
         assert_ulps_eq!(
-            h2.batches[2].posterior(a),
+            h2.batches[2].skills[&a].posterior(),
             Gaussian::new(24.998668, 5.420053),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h2.batches[2].posterior(c),
+            h2.batches[2].skills[&c].posterior(),
             Gaussian::new(25.000532, 5.419827),
             epsilon = 1e-6
         );
@@ -658,17 +658,17 @@ mod tests {
         assert_eq!(h.batches[2].skills[&c].elapsed, 1);
 
         assert_ulps_eq!(
-            h.batches[0].posterior(a),
+            h.batches[0].skills[&a].posterior(),
             Gaussian::new(25.000267, 5.419381),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[0].posterior(b),
+            h.batches[0].skills[&b].posterior(),
             Gaussian::new(24.999465, 5.419425),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[2].posterior(b),
+            h.batches[2].skills[&b].posterior(),
             Gaussian::new(25.000532, 5.419696),
             epsilon = 1e-6
         );
@@ -712,8 +712,8 @@ mod tests {
         );
 
         assert_ulps_eq!(
-            h.batches[0].posterior(b).mu,
-            -1.0 * h.batches[0].posterior(c).mu,
+            h.batches[0].skills[&b].posterior().mu,
+            -1.0 * h.batches[0].skills[&c].posterior().mu,
             epsilon = 1e-6
         );
 
@@ -732,33 +732,33 @@ mod tests {
         assert_ulps_eq!(p_d_m_hat, 0.172432, epsilon = 1e-6);
 
         assert_ulps_eq!(
-            h.batches[0].posterior(a),
-            h.batches[0].posterior(b),
+            h.batches[0].skills[&a].posterior(),
+            h.batches[0].skills[&b].posterior(),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[0].posterior(c),
-            h.batches[0].posterior(d),
+            h.batches[0].skills[&c].posterior(),
+            h.batches[0].skills[&d].posterior(),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[1].posterior(e),
-            h.batches[1].posterior(f),
+            h.batches[1].skills[&e].posterior(),
+            h.batches[1].skills[&f].posterior(),
             epsilon = 1e-6
         );
 
         assert_ulps_eq!(
-            h.batches[0].posterior(a),
+            h.batches[0].skills[&a].posterior(),
             Gaussian::new(4.084902, 5.106919),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[0].posterior(c),
+            h.batches[0].skills[&c].posterior(),
             Gaussian::new(-0.533029, 5.106919),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[2].posterior(e),
+            h.batches[2].skills[&e].posterior(),
             Gaussian::new(-3.551872, 5.154569),
             epsilon = 1e-6
         );
@@ -795,17 +795,17 @@ mod tests {
         assert_eq!(h.batches[2].skills[&c].elapsed, 1);
 
         assert_ulps_eq!(
-            h.batches[0].posterior(a),
+            h.batches[0].skills[&a].posterior(),
             Gaussian::new(0.000000, 1.300610),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[0].posterior(b),
+            h.batches[0].skills[&b].posterior(),
             Gaussian::new(0.000000, 1.300610),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[2].posterior(b),
+            h.batches[2].skills[&b].posterior(),
             Gaussian::new(0.000000, 1.300610),
             epsilon = 1e-6
         );
@@ -832,22 +832,22 @@ mod tests {
         h.convergence(ITERATIONS, EPSILON, false);
 
         assert_ulps_eq!(
-            h.batches[0].posterior(a),
+            h.batches[0].skills[&a].posterior(),
             Gaussian::new(0.000000, 0.931236),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[3].posterior(a),
+            h.batches[3].skills[&a].posterior(),
             Gaussian::new(0.000000, 0.931236),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[3].posterior(b),
+            h.batches[3].skills[&b].posterior(),
             Gaussian::new(0.000000, 0.931236),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[5].posterior(b),
+            h.batches[5].skills[&b].posterior(),
             Gaussian::new(0.000000, 0.931236),
             epsilon = 1e-6
         );
@@ -884,17 +884,17 @@ mod tests {
         assert_eq!(h.batches[2].skills[&c].elapsed, 1);
 
         assert_ulps_eq!(
-            h.batches[0].posterior(a),
+            h.batches[0].skills[&a].posterior(),
             Gaussian::new(0.000000, 1.300610),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[0].posterior(b),
+            h.batches[0].skills[&b].posterior(),
             Gaussian::new(0.000000, 1.300610),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[2].posterior(b),
+            h.batches[2].skills[&b].posterior(),
             Gaussian::new(0.000000, 1.300610),
             epsilon = 1e-6
         );
@@ -921,22 +921,22 @@ mod tests {
         h.convergence(ITERATIONS, EPSILON, false);
 
         assert_ulps_eq!(
-            h.batches[0].posterior(a),
+            h.batches[0].skills[&a].posterior(),
             Gaussian::new(0.000000, 0.931236),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[3].posterior(a),
+            h.batches[3].skills[&a].posterior(),
             Gaussian::new(0.000000, 0.931236),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[3].posterior(b),
+            h.batches[3].skills[&b].posterior(),
             Gaussian::new(0.000000, 0.931236),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
-            h.batches[5].posterior(b),
+            h.batches[5].skills[&b].posterior(),
             Gaussian::new(0.000000, 0.931236),
             epsilon = 1e-6
         );
@@ -1084,20 +1084,20 @@ mod tests {
         h.convergence(ITERATIONS, EPSILON, false);
 
         assert_ulps_eq!(
-            h.batches[0].posterior(b),
-            h.batches[end].posterior(b),
+            h.batches[0].skills[&b].posterior(),
+            h.batches[end].skills[&b].posterior(),
             epsilon = 1e-6
         );
 
         assert_ulps_eq!(
-            h.batches[0].posterior(c),
-            h.batches[end].posterior(c),
+            h.batches[0].skills[&c].posterior(),
+            h.batches[end].skills[&c].posterior(),
             epsilon = 1e-6
         );
 
         assert_ulps_eq!(
-            h.batches[0].posterior(c),
-            h.batches[0].posterior(b),
+            h.batches[0].skills[&c].posterior(),
+            h.batches[0].skills[&b].posterior(),
             epsilon = 1e-6
         );
 
@@ -1172,20 +1172,20 @@ mod tests {
         h.convergence(ITERATIONS, EPSILON, false);
 
         assert_ulps_eq!(
-            h.batches[0].posterior(b),
-            h.batches[end].posterior(b),
+            h.batches[0].skills[&b].posterior(),
+            h.batches[end].skills[&b].posterior(),
             epsilon = 1e-6
         );
 
         assert_ulps_eq!(
-            h.batches[0].posterior(c),
-            h.batches[end].posterior(c),
+            h.batches[0].skills[&c].posterior(),
+            h.batches[end].skills[&c].posterior(),
             epsilon = 1e-6
         );
 
         assert_ulps_eq!(
-            h.batches[0].posterior(c),
-            h.batches[0].posterior(b),
+            h.batches[0].skills[&c].posterior(),
+            h.batches[0].skills[&b].posterior(),
             epsilon = 1e-6
         );
     }
diff --git a/src/lib.rs b/src/lib.rs
index 2f3f1af..dbe72ef 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -4,15 +4,15 @@ use std::collections::HashMap;
 use std::f64::consts::{FRAC_1_SQRT_2, FRAC_2_SQRT_PI, SQRT_2};
 use std::hash::Hash;
 
-mod agent;
+pub mod agent;
 #[cfg(feature = "approx")]
 mod approx;
-mod batch;
+pub mod batch;
 mod game;
-mod gaussian;
+pub mod gaussian;
 mod history;
 mod message;
-mod player;
+pub mod player;
 
 pub use game::Game;
 pub use gaussian::Gaussian;