diff --git a/examples/README.md b/examples/README.md
index 9e38b07..a8882e2 100644
--- a/examples/README.md
+++ b/examples/README.md
@@ -2,6 +2,25 @@
 
 These are some example scripts to demonstrate the various simulations that can be done, and to verify the simulator by reproducing results of already-published works.
 
+Scripts double as tests: each one that has expected output defines a
+`reference_table` (the values its computed `table` is checked against, in the
+same order as the script's columns) and a `tolerance` (the absolute
+comparison tolerance).  `tests/test_examples.py` runs each script and checks
+`table` against `reference_table`.  Where a script reproduces the published
+results, the reference values come from the paper; where it does not (see
+issues #88 and #91), the reference is the script's "Typical result" as a
+regression guard until the discrepancy is fixed.
+
+## Parallel Monte Carlo
+
+Most scripts that run many elections use [joblib](https://joblib.readthedocs.io/) to parallelize embarrassingly parallel simulation loops. The usual pattern is:
+
+- **`batch_size = 100`** — each worker runs this many elections before returning, which cuts down on scheduling overhead.
+- **`Parallel(n_jobs=-3, verbose=5, backend='loky')`** — use all but two CPU cores, print progress, and pin the loky backend so each worker draws elections with its own RNG state (fork-based backends would duplicate the shared RNG across workers).
+- **Aggregation** — workers return partial `Counter` or `defaultdict` results that the main process merges (order does not matter).
+
+These are Monte Carlo verification scripts: they do not fix a global random seed, so exact counts differ from run to run; each worker draws elections independently.
+
 ## Wikipedia
 
 ### Likelihood of a Condorcet cycle
diff --git a/examples/merrill_1984_fig_2c_2d.py b/examples/merrill_1984_fig_2c_2d.py
index 1fb57e0..c514adf 100644
--- a/examples/merrill_1984_fig_2c_2d.py
+++ b/examples/merrill_1984_fig_2c_2d.py
@@ -41,11 +41,11 @@
 running as many simulations, however the Coombs results are consistently
 high.
 """
-import time
 from collections import Counter
 
 import matplotlib.pyplot as plt
 import numpy as np
+from joblib import Parallel, delayed
 from tabulate import tabulate
 
 from elsim.elections import normal_electorate, normed_dist_utilities
@@ -59,6 +59,11 @@
 corr = 0.5
 D = 2
 
+# Simulate more than just one election per worker to improve efficiency
+batch_size = 100
+n_batches = n_elections // batch_size
+assert n_batches * batch_size == n_elections
+
 ranked_methods = {'Plurality': fptp, 'Runoff': runoff, 'Hare': irv,
                   'Borda': borda, 'Coombs': coombs, 'Black': black}
 
@@ -87,14 +92,13 @@
     'Plurality': {2: 100.0, 3: 51.3,  4: 36.2,  5: 21.0,  7: 7.8},
 }
 
-for fig, disp, ymin, orig in (('2.c', 1.0, 50, merrill_fig_2c),
-                              ('2.d', 0.5, 0, merrill_fig_2d)):
+table = {}
 
+def simulate_batch(disp):
+    """Run one batch of elections and return the partial tallies."""
     condorcet_winner_count = {key: Counter() for key in (
         ranked_methods.keys() | rated_methods.keys() | {'CW'})}
-    start_time = time.monotonic()
-
-    for _iteration in range(n_elections):
+    for _iteration in range(batch_size):
         for n_cands in n_cands_list:
             v, c = normal_electorate(n_voters, n_cands, dims=D, corr=corr,
                                      disp=disp)
@@ -115,9 +119,22 @@
                     if method(utilities, tiebreaker='random') == CW:
                         condorcet_winner_count[name][n_cands] += 1
 
-    elapsed_time = time.monotonic() - start_time
-    print('Elapsed:', time.strftime("%H:%M:%S", time.gmtime(elapsed_time)),
-          '\n')
+    return condorcet_winner_count
+
+
+for fig, disp, ymin, orig in (('2.c', 1.0, 50, merrill_fig_2c),
+                              ('2.d', 0.5, 0, merrill_fig_2d)):
+
+    jobs = [delayed(simulate_batch)(disp)] * n_batches
+    print(f'{len(jobs)} tasks total:')
+    results = Parallel(n_jobs=-3, verbose=5, backend='loky')(jobs)
+
+    # Aggregate results
+    condorcet_winner_count = {key: Counter() for key in (
+        ranked_methods.keys() | rated_methods.keys() | {'CW'})}
+    for result in results:
+        for method, counter in result.items():
+            condorcet_winner_count[method].update(counter)
 
     plt.figure(f'Figure {fig}. {n_voters} voters, {n_elections} elections')
     plt.title(f'Figure {fig}: Condorcet Efficiency under Spatial-Model '
@@ -134,7 +151,7 @@
     # Restart color cycle, so result colors match
     plt.gca().set_prop_cycle(None)
 
-    table = []
+    table[fig] = []
 
     # Of those elections with CW, likelihood that method chooses CW
     x_cw, y_cw = zip(*sorted(condorcet_winner_count['CW'].items()))
@@ -143,9 +160,10 @@
         x, y = zip(*sorted(condorcet_winner_count[method].items()))
         CE = np.array(y)/y_cw
         plt.plot(x, CE*100, '-', label=method)
-        table.append([method, *CE*100])
+        table[fig].append([method, *CE*100])
 
-    print(tabulate(table, ["Method", *x], tablefmt="pipe", floatfmt='.1f'))
+    print(tabulate(table[fig], ["Method", *x], tablefmt="pipe",
+                   floatfmt='.1f'))
     print()
 
     plt.plot([], [], 'k:', lw=0.8, label='Merrill')  # Dummy plot for label
@@ -155,3 +173,32 @@
     plt.ylim(ymin, 102)
     plt.xlim(1.8, 7.2)
     plt.show()
+
+# Regression reference from the "Results with 500_000 elections" tables in the
+# docstring, ordered by n_cands_list.  These sims do not reproduce Merrill's
+# published figures exactly (discrepancies up to ~7%, see issue #88), so this
+# is only a regression guard against breaking the current close-enough output
+# until that discrepancy is fixed.
+reference_table = {
+    '2.c': {
+        'Black': (100.0, 100.0, 100.0, 100.0, 100.0, 100.0),
+        'Coombs': (100.0, 99.4, 98.6, 97.8, 96.9, 96.0),
+        'Borda': (100.0, 91.4, 89.2, 87.1, 85.7, 84.6),
+        'Approval': (100.0, 85.9, 79.8, 73.9, 70.1, 66.8),
+        'Hare': (100.0, 94.1, 86.6, 78.9, 71.7, 65.2),
+        'Runoff': (100.0, 94.1, 87.1, 79.7, 72.8, 66.1),
+        'Plurality': (100.0, 80.6, 67.6, 57.4, 49.3, 42.6),
+    },
+    '2.d': {
+        'Black': (100.0, 100.0, 100.0, 100.0, 100.0, 100.0),
+        'Coombs': (100.0, 98.2, 95.9, 93.4, 90.9, 88.4),
+        'Borda': (100.0, 89.2, 86.3, 83.8, 82.1, 80.8),
+        'Approval': (100.0, 84.0, 76.9, 71.5, 67.8, 64.7),
+        'Hare': (100.0, 72.2, 50.3, 35.8, 26.0, 19.7),
+        'Runoff': (100.0, 72.2, 50.6, 35.3, 24.4, 16.9),
+        'Plurality': (100.0, 55.9, 34.7, 21.5, 13.5, 8.5),
+    },
+}
+
+# Absolute tolerance (percentage points) for test_examples.py
+tolerance = 5.0
diff --git a/examples/merrill_1984_fig_2c_2d_updated.py b/examples/merrill_1984_fig_2c_2d_updated.py
index 7a2f480..b8b5f00 100644
--- a/examples/merrill_1984_fig_2c_2d_updated.py
+++ b/examples/merrill_1984_fig_2c_2d_updated.py
@@ -45,11 +45,11 @@
 running as many simulations, however the Coombs results are consistently
 high.
 """
-import time
 from collections import Counter
 
 import matplotlib.pyplot as plt
 import numpy as np
+from joblib import Parallel, delayed
 from tabulate import tabulate
 
 from elsim.elections import normal_electorate, normed_dist_utilities
@@ -64,6 +64,11 @@
 corr = 0.5
 D = 2
 
+# Simulate more than just one election per worker to improve efficiency
+batch_size = 100
+n_batches = n_elections // batch_size
+assert n_batches * batch_size == n_elections
+
 ranked_methods = {'Plurality': fptp, 'Top-2 Runoff': runoff, 'Hare RCV': irv,
                   'Borda': borda, 'Coombs': coombs, 'Condorcet RCV': black}
 
@@ -76,14 +81,13 @@
                      star(honest_normed_scores(utilities, 5), tiebreaker),
                  }
 
-for fig, disp, ymin in (('2.c', 1.0, 50),
-                        ('2.d', 0.5, 0)):
+table = {}
 
+def simulate_batch(disp):
+    """Run one batch of elections and return the partial tallies."""
     condorcet_winner_count = {key: Counter() for key in (
         ranked_methods.keys() | rated_methods.keys() | {'CW'})}
-    start_time = time.monotonic()
-
-    for _iteration in range(n_elections):
+    for _iteration in range(batch_size):
         for n_cands in n_cands_list:
             v, c = normal_electorate(n_voters, n_cands, dims=D, corr=corr,
                                      disp=disp)
@@ -104,16 +108,29 @@
                     if method(utilities, tiebreaker='random') == CW:
                         condorcet_winner_count[name][n_cands] += 1
 
-    elapsed_time = time.monotonic() - start_time
-    print('Elapsed:', time.strftime("%H:%M:%S", time.gmtime(elapsed_time)),
-          '\n')
+    return condorcet_winner_count
+
+
+for fig, disp, ymin in (('2.c', 1.0, 50),
+                        ('2.d', 0.5, 0)):
+
+    jobs = [delayed(simulate_batch)(disp)] * n_batches
+    print(f'{len(jobs)} tasks total:')
+    results = Parallel(n_jobs=-3, verbose=5, backend='loky')(jobs)
+
+    # Aggregate results
+    condorcet_winner_count = {key: Counter() for key in (
+        ranked_methods.keys() | rated_methods.keys() | {'CW'})}
+    for result in results:
+        for method, counter in result.items():
+            condorcet_winner_count[method].update(counter)
 
     plt.figure(f'Figure {fig}. {n_voters} voters, {n_elections} elections',
                figsize=(8, 6.5))
     plt.title(f'Figure {fig}: Condorcet Efficiency under Spatial-Model '
               f'Assumptions [Disp: {disp}]')
 
-    table = []
+    table[fig] = []
 
     # Of those elections with CW, likelihood that method chooses CW
     x_cw, y_cw = zip(*sorted(condorcet_winner_count['CW'].items()))
@@ -122,9 +139,10 @@
         x, y = zip(*sorted(condorcet_winner_count[method].items()))
         CE = np.array(y)/y_cw
         plt.plot(x, CE*100, '-', label=method)
-        table.append([method, *CE*100])
+        table[fig].append([method, *CE*100])
 
-    print(tabulate(table, ["Method", *x], tablefmt="pipe", floatfmt='.1f'))
+    print(tabulate(table[fig], ["Method", *x], tablefmt="pipe",
+                   floatfmt='.1f'))
     print()
 
     plt.legend()
@@ -133,3 +151,36 @@
     plt.ylim(ymin, 102)
     plt.xlim(1.8, 7.2)
     plt.show()
+
+# Regression reference from the "Results with 100_000 elections" tables in the
+# docstring, ordered by n_cands_list.  These sims do not reproduce Merrill's
+# published figures exactly (discrepancies up to ~7%, see issue #88), so this
+# is only a regression guard against breaking the current close-enough output
+# until that discrepancy is fixed.
+reference_table = {
+    '2.c': {
+        'Condorcet RCV': (100.0, 100.0, 100.0, 100.0, 100.0, 100.0),
+        'Coombs': (100.0, 99.4, 98.6, 97.8, 96.9, 96.0),
+        'STAR': (100.0, 97.8, 94.7, 92.1, 89.9, 88.0),
+        'Borda': (100.0, 91.4, 89.2, 87.1, 85.8, 84.8),
+        'Score': (100.0, 88.7, 84.6, 82.7, 81.3, 79.9),
+        'Approval (opt.)': (100.0, 86.0, 79.7, 73.9, 70.5, 67.0),
+        'Hare RCV': (100.0, 94.1, 86.6, 79.0, 71.3, 65.1),
+        'Top-2 Runoff': (100.0, 94.1, 87.1, 79.9, 72.8, 65.9),
+        'Plurality': (100.0, 80.6, 67.8, 57.3, 49.1, 42.7),
+    },
+    '2.d': {
+        'Condorcet RCV': (100.0, 100.0, 100.0, 100.0, 100.0, 100.0),
+        'Coombs': (100.0, 98.3, 95.7, 93.4, 90.8, 88.5),
+        'STAR': (100.0, 96.5, 91.5, 87.2, 83.1, 79.8),
+        'Borda': (100.0, 89.2, 86.4, 83.9, 82.2, 80.7),
+        'Score': (100.0, 86.2, 80.5, 77.4, 74.5, 72.1),
+        'Approval (opt.)': (100.0, 83.7, 76.9, 71.6, 67.6, 64.6),
+        'Hare RCV': (100.0, 72.4, 50.4, 35.6, 26.2, 19.8),
+        'Top-2 Runoff': (100.0, 72.4, 50.5, 35.3, 24.5, 17.1),
+        'Plurality': (100.0, 56.3, 34.8, 21.5, 13.5, 8.6),
+    },
+}
+
+# Absolute tolerance (percentage points) for test_examples.py
+tolerance = 5.0
diff --git a/examples/merrill_1984_fig_4a_4b.py b/examples/merrill_1984_fig_4a_4b.py
index 2a75ce8..ea8279d 100644
--- a/examples/merrill_1984_fig_4a_4b.py
+++ b/examples/merrill_1984_fig_4a_4b.py
@@ -40,12 +40,12 @@
 discrepancies.  It is smoother, so maybe the original just had lower number of
 simulations.
 """
-import time
 from collections import Counter
 from random import randint
 
 import matplotlib.pyplot as plt
 import numpy as np
+from joblib import Parallel, delayed
 from tabulate import tabulate
 
 from elsim.elections import normal_electorate, normed_dist_utilities
@@ -59,6 +59,11 @@
 corr = 0.5
 D = 2
 
+# Simulate more than just one election per worker to improve efficiency
+batch_size = 100
+n_batches = n_elections // batch_size
+assert n_batches * batch_size == n_elections
+
 ranked_methods = {'Plurality': fptp, 'Runoff': runoff, 'Hare': irv,
                   'Borda': borda, 'Coombs': coombs, 'Black': black}
 
@@ -87,15 +92,14 @@
     'Plurality': {2: 100.0, 3: 41.1, 4: 27.0, 5: -1.0, 7: -9},
 }
 
-for fig, disp, ymin, orig in (('4.a', 1.0, 55, merrill_fig_4a),
-                              ('4.b', 0.5, 0, merrill_fig_4b)):
+table = {}
 
+def simulate_batch(disp):
+    """Run one batch of elections and return the partial tallies."""
     utility_sums = {key: Counter() for key in (ranked_methods.keys() |
                                                rated_methods.keys() |
                                                {'SU max', 'RW'})}
-    start_time = time.monotonic()
-
-    for _iteration in range(n_elections):
+    for _iteration in range(batch_size):
         for n_cands in n_cands_list:
             v, c = normal_electorate(n_voters, n_cands, dims=D, corr=corr,
                                      disp=disp)
@@ -114,9 +118,23 @@
                 winner = method(rankings, tiebreaker='random')
                 utility_sums[name][n_cands] += utilities.sum(axis=0)[winner]
 
-    elapsed_time = time.monotonic() - start_time
-    print('Elapsed:', time.strftime("%H:%M:%S", time.gmtime(elapsed_time)),
-          '\n')
+    return utility_sums
+
+
+for fig, disp, ymin, orig in (('4.a', 1.0, 55, merrill_fig_4a),
+                              ('4.b', 0.5, 0, merrill_fig_4b)):
+
+    jobs = [delayed(simulate_batch)(disp)] * n_batches
+    print(f'{len(jobs)} tasks total:')
+    results = Parallel(n_jobs=-3, verbose=5, backend='loky')(jobs)
+
+    # Aggregate results
+    utility_sums = {key: Counter() for key in (ranked_methods.keys() |
+                                               rated_methods.keys() |
+                                               {'SU max', 'RW'})}
+    for result in results:
+        for method, counter in result.items():
+            utility_sums[method].update(counter)
 
     plt.figure(f'Figure {fig}. {n_voters} voters, {n_elections} elections')
     plt.title(f'Figure {fig}: Social Utility Efficiency under Spatial-Model '
@@ -133,7 +151,7 @@
     # Restart color cycle, so result colors match
     plt.gca().set_prop_cycle(None)
 
-    table = []
+    table[fig] = []
 
     # Calculate Social Utility Efficiency from summed utilities
     x_uw, y_uw = zip(*sorted(utility_sums['SU max'].items()))
@@ -143,9 +161,10 @@
         x, y = zip(*sorted(utility_sums[method].items()))
         SUE = (np.array(y) - y_rw) / (np.array(y_uw) - y_rw)
         plt.plot(x, SUE*100, '-', label=method)
-        table.append([method, *SUE*100])
+        table[fig].append([method, *SUE*100])
 
-    print(tabulate(table, ["Method", *x], tablefmt="pipe", floatfmt='.1f'))
+    print(tabulate(table[fig], ["Method", *x], tablefmt="pipe",
+                   floatfmt='.1f'))
     print()
 
     plt.plot([], [], 'k:', lw=0.8, label='Merrill')  # Dummy plot for label
@@ -155,3 +174,32 @@
     plt.ylim(ymin, 102)
     plt.xlim(1.8, 7.2)
     plt.show()
+
+# Regression reference from the "Results with 500_000 elections" tables in the
+# docstring, ordered by n_cands_list.  These sims do not reproduce Merrill's
+# published figures exactly (discrepancies up to ~7%, see issue #88), so this
+# is only a regression guard against breaking the current close-enough output
+# until that discrepancy is fixed.
+reference_table = {
+    '4.a': {
+        'Black': (100.0, 97.2, 97.1, 97.3, 97.6, 97.8),
+        'Coombs': (100.0, 97.1, 96.8, 97.0, 97.2, 97.4),
+        'Borda': (100.0, 98.7, 98.2, 97.9, 97.7, 97.6),
+        'Approval': (100.0, 98.7, 97.3, 96.2, 95.6, 95.2),
+        'Hare': (100.0, 94.2, 92.6, 91.7, 91.0, 90.3),
+        'Runoff': (100.0, 94.2, 92.0, 90.4, 88.9, 87.4),
+        'Plurality': (100.0, 84.7, 77.1, 72.1, 68.1, 64.8),
+    },
+    '4.b': {
+        'Black': (100.0, 95.5, 95.2, 95.5, 95.8, 96.2),
+        'Coombs': (100.0, 94.9, 94.1, 94.0, 94.0, 94.1),
+        'Borda': (100.0, 97.9, 97.1, 96.6, 96.4, 96.3),
+        'Approval': (100.0, 98.6, 96.7, 95.6, 94.9, 94.5),
+        'Hare': (100.0, 70.2, 55.9, 46.7, 39.7, 34.6),
+        'Runoff': (100.0, 70.2, 51.7, 36.9, 24.3, 13.5),
+        'Plurality': (100.0, 50.1, 23.7, 4.3, -11.8, -25.1),
+    },
+}
+
+# Absolute tolerance (percentage points) for test_examples.py
+tolerance = 5.0
diff --git a/examples/merrill_1984_fig_4a_4b_updated.py b/examples/merrill_1984_fig_4a_4b_updated.py
index 28feb39..bad51c1 100644
--- a/examples/merrill_1984_fig_4a_4b_updated.py
+++ b/examples/merrill_1984_fig_4a_4b_updated.py
@@ -44,12 +44,12 @@
 discrepancies.  It is smoother, so maybe the original just had lower number of
 simulations.
 """
-import time
 from collections import Counter
 from random import randint
 
 import matplotlib.pyplot as plt
 import numpy as np
+from joblib import Parallel, delayed
 from tabulate import tabulate
 
 from elsim.elections import normal_electorate, normed_dist_utilities
@@ -64,6 +64,11 @@
 corr = 0.5
 D = 2
 
+# Simulate more than just one election per worker to improve efficiency
+batch_size = 100
+n_batches = n_elections // batch_size
+assert n_batches * batch_size == n_elections
+
 ranked_methods = {'Plurality': fptp, 'Top-2 Runoff': runoff, 'Hare RCV': irv,
                   'Borda': borda, 'Coombs': coombs, 'Condorcet RCV': black}
 
@@ -78,15 +83,14 @@
                           tiebreaker),
                  }
 
-for fig, disp, _ymin in (('4.a', 1.0, 55),
-                        ('4.b', 0.5, 0)):
+table = {}
 
+def simulate_batch(disp):
+    """Run one batch of elections and return the partial tallies."""
     utility_sums = {key: Counter() for key in (ranked_methods.keys() |
                                                rated_methods.keys() |
                                                {'SU max', 'RW'})}
-    start_time = time.monotonic()
-
-    for _iteration in range(n_elections):
+    for _iteration in range(batch_size):
         for n_cands in n_cands_list:
             v, c = normal_electorate(n_voters, n_cands, dims=D, corr=corr,
                                      disp=disp)
@@ -105,16 +109,30 @@
                 winner = method(rankings, tiebreaker='random')
                 utility_sums[name][n_cands] += utilities.sum(axis=0)[winner]
 
-    elapsed_time = time.monotonic() - start_time
-    print('Elapsed:', time.strftime("%H:%M:%S", time.gmtime(elapsed_time)),
-          '\n')
+    return utility_sums
+
+
+for fig, disp, _ymin in (('4.a', 1.0, 55),
+                        ('4.b', 0.5, 0)):
+
+    jobs = [delayed(simulate_batch)(disp)] * n_batches
+    print(f'{len(jobs)} tasks total:')
+    results = Parallel(n_jobs=-3, verbose=5, backend='loky')(jobs)
+
+    # Aggregate results
+    utility_sums = {key: Counter() for key in (ranked_methods.keys() |
+                                               rated_methods.keys() |
+                                               {'SU max', 'RW'})}
+    for result in results:
+        for method, counter in result.items():
+            utility_sums[method].update(counter)
 
     plt.figure(f'Figure {fig}. {n_voters} voters, {n_elections} elections',
                figsize=(8, 6.5))
     plt.title(f'Figure {fig}: Social Utility Efficiency under Spatial-Model '
               f'Assumptions [Disp: {disp}]')
 
-    table = []
+    table[fig] = []
 
     # Calculate Social Utility Efficiency from summed utilities
     x_uw, y_uw = zip(*sorted(utility_sums['SU max'].items()))
@@ -124,9 +142,10 @@
         x, y = zip(*sorted(utility_sums[method].items()))
         SUE = (np.array(y) - y_rw) / (np.array(y_uw) - y_rw)
         plt.plot(x, SUE*100, '-', label=method)
-        table.append([method, *SUE*100])
+        table[fig].append([method, *SUE*100])
 
-    print(tabulate(table, ["Method", *x], tablefmt="pipe", floatfmt='.1f'))
+    print(tabulate(table[fig], ["Method", *x], tablefmt="pipe",
+                   floatfmt='.1f'))
     print()
 
     plt.legend()
@@ -135,3 +154,36 @@
     plt.ylim(85, 100.5)  # or ymin
     plt.xlim(1.8, 7.2)
     plt.show()
+
+# Regression reference from the "Results with 100_000 elections" tables in the
+# docstring, ordered by n_cands_list.  These sims do not reproduce Merrill's
+# published figures exactly (discrepancies up to ~7%, see issue #88), so this
+# is only a regression guard against breaking the current close-enough output
+# until that discrepancy is fixed.
+reference_table = {
+    '4.a': {
+        'Score': (100.0, 100.0, 99.9, 99.9, 99.9, 99.8),
+        'STAR': (100.0, 97.5, 97.8, 98.3, 98.6, 98.8),
+        'Borda': (100.0, 98.8, 98.2, 97.9, 97.6, 97.6),
+        'Condorcet RCV': (100.0, 97.1, 97.0, 97.3, 97.6, 97.8),
+        'Coombs': (100.0, 97.0, 96.8, 97.0, 97.2, 97.4),
+        'Approval (opt.)': (100.0, 98.7, 97.3, 96.2, 95.5, 95.2),
+        'Hare RCV': (100.0, 94.2, 92.5, 91.6, 91.0, 90.4),
+        'Top-2 Runoff': (100.0, 94.2, 91.9, 90.4, 88.9, 87.5),
+        'Plurality': (100.0, 84.8, 77.4, 72.0, 68.2, 64.9),
+    },
+    '4.b': {
+        'Score': (100.0, 100.0, 99.9, 99.7, 99.5, 99.3),
+        'STAR': (100.0, 96.2, 96.7, 97.2, 97.6, 97.8),
+        'Borda': (100.0, 97.9, 97.1, 96.6, 96.4, 96.3),
+        'Condorcet RCV': (100.0, 95.5, 95.2, 95.5, 95.8, 96.2),
+        'Coombs': (100.0, 95.0, 94.2, 94.1, 94.0, 94.1),
+        'Approval (opt.)': (100.0, 98.6, 96.7, 95.5, 94.9, 94.6),
+        'Hare RCV': (100.0, 70.4, 55.9, 46.7, 39.5, 35.0),
+        'Top-2 Runoff': (100.0, 70.4, 51.8, 37.3, 23.9, 13.6),
+        'Plurality': (100.0, 50.5, 23.9, 4.7, -12.2, -24.7),
+    },
+}
+
+# Absolute tolerance (percentage points) for test_examples.py
+tolerance = 5.0
diff --git a/examples/merrill_1984_table_1_fig_1.py b/examples/merrill_1984_table_1_fig_1.py
index 951b0ea..ab740de 100644
--- a/examples/merrill_1984_table_1_fig_1.py
+++ b/examples/merrill_1984_table_1_fig_1.py
@@ -24,11 +24,11 @@
 | SU max    | 100.0 |  84.1 |  79.6 |  78.4 |  77.3 |  77.5 |
 | CW        | 100.0 |  91.7 |  83.1 |  75.6 |  64.3 |  52.9 |
 """
-import time
 from collections import Counter
 
 import matplotlib.pyplot as plt
 import numpy as np
+from joblib import Parallel, delayed
 from tabulate import tabulate
 
 from elsim.elections import random_utilities
@@ -40,6 +40,11 @@
 n_voters = 25
 n_cands_list = (2, 3, 4, 5, 7, 10)
 
+# Simulate more than just one election per worker to improve efficiency
+batch_size = 100
+n_batches = n_elections // batch_size
+assert n_batches * batch_size == n_elections
+
 ranked_methods = {'Plurality': fptp, 'Runoff': runoff, 'Hare': irv,
                   'Borda': borda, 'Coombs': coombs, 'Black': black}
 
@@ -47,45 +52,62 @@
                  'Approval': lambda utilities, tiebreaker:
                      approval(approval_optimal(utilities), tiebreaker)}
 
-condorcet_winner_count = {key: Counter() for key in (
-    ranked_methods.keys() | rated_methods.keys() | {'CW'})}
 
-start_time = time.monotonic()
+def simulate_batch():
+    """Run one batch of elections and return the partial tallies."""
+    condorcet_winner_count = {key: Counter() for key in (
+        ranked_methods.keys() | rated_methods.keys() | {'CW'})}
+
+    for _iteration in range(batch_size):
+        for n_cands in n_cands_list:
+            utilities = random_utilities(n_voters, n_cands)
+
+            """
+            "Simulated utilities were normalized by range, that is, each
+            voter's set of utilities were linearly expanded so that the highest
+            and lowest utilities for each voter were 1 and 0, respectively."
+
+            This is necessary for the SU Maximizer results to match Merrill's.
+            """
+            utilities -= utilities.min(1)[:, np.newaxis]
+            utilities /= utilities.max(1)[:, np.newaxis]
 
-for _iteration in range(n_elections):
-    for n_cands in n_cands_list:
-        utilities = random_utilities(n_voters, n_cands)
+            rankings = honest_rankings(utilities)
 
-        """
-        "Simulated utilities were normalized by range, that is, each voter's
-        set of utilities were linearly expanded so that the highest and lowest
-        utilities for each voter were 1 and 0, respectively."
+            # If there is a Condorcet winner, analyze election, otherwise skip
+            # it
+            CW = condorcet(rankings)
+            if CW is not None:
+                condorcet_winner_count['CW'][n_cands] += 1
 
-        This is necessary for the SU Maximizer results to match Merrill's.
-        """
-        utilities -= utilities.min(1)[:, np.newaxis]
-        utilities /= utilities.max(1)[:, np.newaxis]
+                for name, method in ranked_methods.items():
+                    if method(rankings, tiebreaker='random') == CW:
+                        condorcet_winner_count[name][n_cands] += 1
 
-        rankings = honest_rankings(utilities)
+                for name, method in rated_methods.items():
+                    if method(utilities, tiebreaker='random') == CW:
+                        condorcet_winner_count[name][n_cands] += 1
 
-        # If there is a Condorcet winner, analyze election, otherwise skip it
-        CW = condorcet(rankings)
-        if CW is not None:
-            condorcet_winner_count['CW'][n_cands] += 1
+    return condorcet_winner_count
 
-            for name, method in ranked_methods.items():
-                if method(rankings, tiebreaker='random') == CW:
-                    condorcet_winner_count[name][n_cands] += 1
 
-            for name, method in rated_methods.items():
-                if method(utilities, tiebreaker='random') == CW:
-                    condorcet_winner_count[name][n_cands] += 1
 
-elapsed_time = time.monotonic() - start_time
-print('Elapsed:', time.strftime("%H:%M:%S", time.gmtime(elapsed_time)), '\n')
+jobs = [delayed(simulate_batch)()] * n_batches
+print(f'{len(jobs)} tasks total:')
+results = Parallel(n_jobs=-3, verbose=5, backend='loky')(jobs)
 
-# Plot Merrill's results as dotted lines for comparison
-merrill_table_1 = {
+# Aggregate results
+condorcet_winner_count = {key: Counter() for key in (
+    ranked_methods.keys() | rated_methods.keys() | {'CW'})}
+for result in results:
+    for method, counter in result.items():
+        condorcet_winner_count[method].update(counter)
+
+
+# Reference values from Merrill's published Table 1, used to plot his results
+# as dotted lines for comparison.  This script reproduces them within ~2 pp,
+# so test_examples.py checks the computed table against these.
+reference_table = {
     'Plurality': {2: 100.0, 3:  79.1, 4:  69.4, 5:  62.1, 7:  52.0, 10:  42.6},
     'Runoff':    {2: 100.0, 3:  96.2, 4:  90.1, 5:  83.6, 7:  73.5, 10:  61.3},
     'Hare':      {2: 100.0, 3:  96.2, 4:  92.7, 5:  89.1, 7:  84.8, 10:  77.9},
@@ -97,11 +119,14 @@
     'CW':        {2: 100.0, 3:  91.6, 4:  83.4, 5:  75.8, 7:  64.3, 10:  52.5},
 }
 
+# Absolute tolerance (percentage points) for test_examples.py
+tolerance = 4.0
+
 plt.figure(f'Figure 1. {n_voters} voters, {n_elections} elections')
 plt.title('Figure 1: Condorcet Efficiencies for a Random Society')
 for method in ('Plurality', 'Runoff', 'Hare', 'Approval', 'Borda', 'Coombs',
                'Black'):
-    x, y = zip(*sorted(merrill_table_1[method].items()))
+    x, y = zip(*sorted(reference_table[method].items()))
     plt.plot(x, y, ':', lw=0.8)
 
 # Restart color cycle, so result colors match
diff --git a/examples/merrill_1984_table_2.py b/examples/merrill_1984_table_2.py
index 2948f6a..a149cd9 100644
--- a/examples/merrill_1984_table_2.py
+++ b/examples/merrill_1984_table_2.py
@@ -31,10 +31,10 @@
 Many of these values match the paper closely, but some are consistently off by
 up to 4%.
 """
-import time
 from collections import Counter
 
 import numpy as np
+from joblib import Parallel, delayed
 from tabulate import tabulate
 
 from elsim.elections import normal_electorate, normed_dist_utilities
@@ -46,6 +46,11 @@
 n_voters = 201
 n_cands = 5
 
+# Simulate more than just one election per worker to improve efficiency
+batch_size = 100
+n_batches = n_elections // batch_size
+assert n_batches * batch_size == n_elections
+
 ranked_methods = {'Plurality': fptp, 'Runoff': runoff, 'Hare': irv,
                   'Borda': borda, 'Coombs': coombs, 'Black': black}
 
@@ -53,8 +58,6 @@
                  'Approval': lambda utilities, tiebreaker:
                      approval(approval_optimal(utilities), tiebreaker)}
 
-start_time = time.monotonic()
-
 #             disp, corr, D
 conditions = ((1.0, 0.5, 2),
               (1.0, 0.5, 4),
@@ -66,14 +69,11 @@
               (0.5, 0.0, 4),
               )
 
-results = []
-
-for disp, corr, D in conditions:
-    print(disp, corr, D)
 
+def simulate_batch(disp, corr, D):
+    """Run one batch of elections and return the partial tallies."""
     condorcet_winner_count = Counter()
-
-    for _iteration in range(n_elections):
+    for _iteration in range(batch_size):
         v, c = normal_electorate(n_voters, n_cands, dims=D, corr=corr,
                                  disp=disp)
 
@@ -103,10 +103,26 @@
                 if method(utilities, tiebreaker='random') == CW:
                     condorcet_winner_count[name] += 1
 
+    return condorcet_winner_count
+
+
+
+results = []
+
+for disp, corr, D in conditions:
+    print(disp, corr, D)
+
+    jobs = [delayed(simulate_batch)(disp, corr, D)] * n_batches
+    print(f'{len(jobs)} tasks total:')
+    batch_results = Parallel(n_jobs=-3, verbose=5, backend='loky')(jobs)
+
+    # Aggregate results for this condition
+    condorcet_winner_count = Counter()
+    for result in batch_results:
+        condorcet_winner_count.update(result)
+
     results.append(condorcet_winner_count)
 
-elapsed_time = time.monotonic() - start_time
-print('Elapsed:', time.strftime("%H:%M:%S", time.gmtime(elapsed_time)), '\n')
 
 # Neither Tabulate nor Markdown support column span or multiple headers, but
 # at least this prints to plain text in a readable way.
@@ -124,3 +140,23 @@
 table.append(['CW', *(y_cw / n_elections * 100)])
 
 print(tabulate(table, header, tablefmt="pipe", floatfmt='.1f'))
+
+# Regression reference from the "Typical result" table in the docstring,
+# ordered by `conditions`.  The published Merrill (1984) Table 2 is not
+# reproduced (some values are off by up to ~5 pp; see issue #88), so this is
+# only a regression guard against breaking the current close-enough output
+# until that discrepancy is fixed.
+reference_table = {
+    'Plurality': (57.5, 65.8, 62.2, 78.4, 21.7, 24.4, 27.2, 41.3),
+    'Runoff': (80.1, 87.3, 81.6, 93.6, 35.4, 42.2, 41.5, 61.5),
+    'Hare': (79.2, 86.7, 84.0, 95.4, 35.9, 46.8, 41.0, 69.9),
+    'Approval': (73.8, 77.8, 76.9, 85.4, 71.5, 76.4, 73.8, 82.7),
+    'Borda': (87.1, 89.3, 88.2, 92.3, 83.7, 86.3, 85.2, 89.4),
+    'Coombs': (97.8, 97.3, 97.9, 98.2, 93.5, 92.3, 93.8, 94.5),
+    'Black': (100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0),
+    'SU max': (82.9, 85.8, 85.3, 90.8, 78.1, 81.5, 80.8, 87.1),
+    'CW': (99.7, 99.7, 99.7, 99.6, 98.9, 98.6, 98.7, 98.5),
+}
+
+# Absolute tolerance (percentage points) for test_examples.py
+tolerance = 3.5
diff --git a/examples/merrill_1984_table_3_fig_3.py b/examples/merrill_1984_table_3_fig_3.py
index d58804c..807edb7 100644
--- a/examples/merrill_1984_table_3_fig_3.py
+++ b/examples/merrill_1984_table_3_fig_3.py
@@ -22,11 +22,11 @@
 | Coombs    | 100.0 | 90.2 | 86.8 | 85.2 | 84.0 | 82.9 |
 | Black     | 100.0 | 92.9 | 92.0 | 92.1 | 93.2 | 94.6 |
 """
-import time
 from collections import Counter
 
 import matplotlib.pyplot as plt
 import numpy as np
+from joblib import Parallel, delayed
 from tabulate import tabulate
 
 from elsim.elections import random_utilities
@@ -38,49 +38,69 @@
 n_voters = 25
 n_cands_list = (2, 3, 4, 5, 7, 10)
 
+# Simulate more than just one election per worker to improve efficiency
+batch_size = 100
+n_batches = n_elections // batch_size
+assert n_batches * batch_size == n_elections
+
 ranked_methods = {'Plurality': fptp, 'Runoff': runoff, 'Hare': irv,
                   'Borda': borda, 'Coombs': coombs, 'Black': black}
 
 rated_methods = {'Approval': lambda utilities, tiebreaker:
                  approval(approval_optimal(utilities), tiebreaker)}
 
-utility_sums = {key: Counter() for key in (ranked_methods.keys() |
-                                           rated_methods.keys() | {'UW'})}
 
-start_time = time.monotonic()
+def simulate_batch():
+    """Run one batch of elections and return the partial tallies."""
+    utility_sums = {key: Counter() for key in (ranked_methods.keys() |
+                                               rated_methods.keys() | {'UW'})}
+
+    for _iteration in range(batch_size):
+        for n_cands in n_cands_list:
+            utilities = random_utilities(n_voters, n_cands)
+
+            """
+            "Simulated utilities were normalized by range, that is, each
+            voter's set of utilities were linearly expanded so that the highest
+            and lowest utilities for each voter were 1 and 0, respectively."
+            """
+            # TODO: Try the Standard Score normalization too?
+            utilities -= utilities.min(1)[:, np.newaxis]
+            utilities /= utilities.max(1)[:, np.newaxis]
 
-for _iteration in range(n_elections):
-    for n_cands in n_cands_list:
-        utilities = random_utilities(n_voters, n_cands)
+            # Find the social utility winner and accumulate utilities
+            UW = utility_winner(utilities)
+            utility_sums['UW'][n_cands] += utilities.sum(axis=0)[UW]
 
-        """
-        "Simulated utilities were normalized by range, that is, each voter's
-        set of utilities were linearly expanded so that the highest and lowest
-        utilities for each voter were 1 and 0, respectively."
-        """
-        # TODO: Try the Standard Score normalization too?
-        utilities -= utilities.min(1)[:, np.newaxis]
-        utilities /= utilities.max(1)[:, np.newaxis]
+            for name, method in rated_methods.items():
+                winner = method(utilities, tiebreaker='random')
+                utility_sums[name][n_cands] += utilities.sum(axis=0)[winner]
 
-        # Find the social utility winner and accumulate utilities
-        UW = utility_winner(utilities)
-        utility_sums['UW'][n_cands] += utilities.sum(axis=0)[UW]
+            rankings = honest_rankings(utilities)
+            for name, method in ranked_methods.items():
+                winner = method(rankings, tiebreaker='random')
+                utility_sums[name][n_cands] += utilities.sum(axis=0)[winner]
 
-        for name, method in rated_methods.items():
-            winner = method(utilities, tiebreaker='random')
-            utility_sums[name][n_cands] += utilities.sum(axis=0)[winner]
+    return utility_sums
 
-        rankings = honest_rankings(utilities)
-        for name, method in ranked_methods.items():
-            winner = method(rankings, tiebreaker='random')
-            utility_sums[name][n_cands] += utilities.sum(axis=0)[winner]
 
 
-elapsed_time = time.monotonic() - start_time
-print('Elapsed:', time.strftime("%H:%M:%S", time.gmtime(elapsed_time)), '\n')
+jobs = [delayed(simulate_batch)()] * n_batches
+print(f'{len(jobs)} tasks total:')
+results = Parallel(n_jobs=-3, verbose=5, backend='loky')(jobs)
 
-# Plot Merrill's results as dotted lines for comparison
-merrill_table_1 = {
+# Aggregate results
+utility_sums = {key: Counter() for key in (ranked_methods.keys() |
+                                           rated_methods.keys() | {'UW'})}
+for result in results:
+    for method, counter in result.items():
+        utility_sums[method].update(counter)
+
+
+# Reference values from Merrill's published Table 3, used to plot his results
+# as dotted lines for comparison.  This script reproduces them within ~2 pp,
+# so test_examples.py checks the computed table against these.
+reference_table = {
     'Plurality': {2: 100.0, 3: 83.0, 4: 75.0, 5: 69.2, 7: 62.8, 10: 53.3},
     'Runoff':    {2: 100.0, 3: 89.5, 4: 83.8, 5: 80.5, 7: 75.6, 10: 67.6},
     'Hare':      {2: 100.0, 3: 89.5, 4: 84.7, 5: 82.4, 7: 80.5, 10: 74.9},
@@ -90,11 +110,14 @@
     'Black':     {2: 100.0, 3: 93.1, 4: 91.9, 5: 92.0, 7: 93.1, 10: 94.3},
 }
 
+# Absolute tolerance (percentage points) for test_examples.py
+tolerance = 3.5
+
 plt.figure(f'Figure 3. {n_voters} voters, {n_elections} elections')
 plt.title('Figure 3: Efficiencies for Social Utility for a Random Society')
 for method in ('Plurality', 'Runoff', 'Hare', 'Approval', 'Borda', 'Coombs',
                'Black'):
-    x, y = zip(*sorted(merrill_table_1[method].items()))
+    x, y = zip(*sorted(reference_table[method].items()))
     plt.plot(x, y, ':', lw=0.8)
 
 # Restart color cycle, so result colors match
diff --git a/examples/merrill_1984_table_4.py b/examples/merrill_1984_table_4.py
index 721c8ff..a64627f 100644
--- a/examples/merrill_1984_table_4.py
+++ b/examples/merrill_1984_table_4.py
@@ -29,11 +29,11 @@
 Many of these values match the paper closely, but some are consistently off by
 up to 9%.
 """
-import time
 from collections import Counter
 from random import randint
 
 import numpy as np
+from joblib import Parallel, delayed
 from tabulate import tabulate
 
 from elsim.elections import normal_electorate, normed_dist_utilities
@@ -45,6 +45,11 @@
 n_voters = 201
 n_cands = 5
 
+# Simulate more than just one election per worker to improve efficiency
+batch_size = 100
+n_batches = n_elections // batch_size
+assert n_batches * batch_size == n_elections
+
 ranked_methods = {'Plurality': fptp, 'Runoff': runoff, 'Hare': irv,
                   'Borda': borda, 'Coombs': coombs, 'Black': black}
 
@@ -52,8 +57,6 @@
                  'Approval': lambda utilities, tiebreaker:
                      approval(approval_optimal(utilities), tiebreaker)}
 
-start_time = time.monotonic()
-
 #             disp, corr, D
 conditions = ((1.0, 0.5, 2),
               (1.0, 0.5, 4),
@@ -65,14 +68,11 @@
               (0.5, 0.0, 4),
               )
 
-results = []
-
-for disp, corr, D in conditions:
-    print(disp, corr, D)
 
+def simulate_batch(disp, corr, D):
+    """Run one batch of elections and return the partial tallies."""
     utility_sums = Counter()
-
-    for _iteration in range(n_elections):
+    for _iteration in range(batch_size):
         v, c = normal_electorate(n_voters, n_cands, dims=D, corr=corr,
                                  disp=disp)
 
@@ -101,10 +101,26 @@
             winner = method(rankings, tiebreaker='random')
             utility_sums[name] += utilities.sum(axis=0)[winner]
 
+    return utility_sums
+
+
+
+results = []
+
+for disp, corr, D in conditions:
+    print(disp, corr, D)
+
+    jobs = [delayed(simulate_batch)(disp, corr, D)] * n_batches
+    print(f'{len(jobs)} tasks total:')
+    batch_results = Parallel(n_jobs=-3, verbose=5, backend='loky')(jobs)
+
+    # Aggregate results for this condition
+    utility_sums = Counter()
+    for result in batch_results:
+        utility_sums.update(result)
+
     results.append(utility_sums)
 
-elapsed_time = time.monotonic() - start_time
-print('Elapsed:', time.strftime("%H:%M:%S", time.gmtime(elapsed_time)), '\n')
 
 # Neither Tabulate nor Markdown support column span or multiple headers, but
 # at least this prints to plain text in a readable way.
@@ -121,3 +137,21 @@
     SUE = (y - y_rw)/(y_uw - y_rw)
     table.append([method, *(SUE*100)])
 print(tabulate(table, header, tablefmt="pipe", floatfmt='.1f'))
+
+# Regression reference from the "Typical result" table in the docstring,
+# ordered by `conditions`.  The published Merrill (1984) Table 4 is not
+# reproduced (some values are off by up to ~9 pp; see issue #88), so this is
+# only a regression guard against breaking the current close-enough output
+# until that discrepancy is fixed.
+reference_table = {
+    'Plurality': (72.1, 79.1, 80.4, 92.4, 4.0, 6.3, 25.2, 52.9),
+    'Runoff': (90.5, 94.2, 92.0, 97.5, 36.6, 43.6, 53.3, 75.3),
+    'Hare': (91.7, 94.7, 94.3, 98.4, 46.4, 57.7, 58.7, 83.6),
+    'Approval': (96.2, 97.0, 96.8, 98.5, 95.6, 96.8, 95.8, 98.0),
+    'Borda': (97.8, 98.6, 98.3, 99.4, 96.6, 97.7, 97.4, 99.0),
+    'Coombs': (97.0, 97.5, 97.7, 98.7, 94.0, 94.3, 95.0, 96.7),
+    'Black': (97.3, 97.8, 98.0, 99.0, 95.5, 96.1, 96.5, 98.0),
+}
+
+# Absolute tolerance (percentage points) for test_examples.py
+tolerance = 3.5
diff --git a/examples/weber_1977_effectiveness_table.py b/examples/weber_1977_effectiveness_table.py
index 3fcf53f..7e5e48f 100644
--- a/examples/weber_1977_effectiveness_table.py
+++ b/examples/weber_1977_effectiveness_table.py
@@ -21,11 +21,11 @@
 """
 # TODO: Standard is consistently ~1% high, while Borda is very accurate
 # TODO: Best Vote-for-or-against-k is not implemented yet
-import time
 from collections import Counter
 
 import matplotlib.pyplot as plt
 import numpy as np
+from joblib import Parallel, delayed
 from tabulate import tabulate
 
 from elsim.elections import random_utilities
@@ -37,35 +37,54 @@
 n_voters = 1_000
 n_cands_list = (2, 3, 4, 5, 6, 10, 255)
 
+# Simulate more than just one election per worker to improve efficiency
+batch_size = 100
+n_batches = n_elections // batch_size
+assert n_batches * batch_size == n_elections
+
 ranked_methods = {'Standard': fptp, 'Borda': borda}
 
 rated_methods = {'Vote-for-half': lambda utilities, tiebreaker:
                  approval(vote_for_k(utilities, 'half'), tiebreaker)}
 
-utility_sums = {key: Counter() for key in (ranked_methods.keys() |
-                                           rated_methods.keys() | {'UW'})}
 
-start_time = time.monotonic()
+def simulate_batch():
+    """Run one batch of elections and return the partial tallies."""
+    utility_sums = {key: Counter() for key in (ranked_methods.keys() |
+                                               rated_methods.keys() | {'UW'})}
+
+    for _iteration in range(batch_size):
+        for n_cands in n_cands_list:
+            utilities = random_utilities(n_voters, n_cands)
+
+            # Find the social utility winner and accumulate utilities
+            UW = utility_winner(utilities)
+            utility_sums['UW'][n_cands] += utilities.sum(axis=0)[UW]
+
+            for name, method in rated_methods.items():
+                winner = method(utilities, tiebreaker='random')
+                utility_sums[name][n_cands] += utilities.sum(axis=0)[winner]
 
-for _iteration in range(n_elections):
-    for n_cands in n_cands_list:
-        utilities = random_utilities(n_voters, n_cands)
+            rankings = honest_rankings(utilities)
+            for name, method in ranked_methods.items():
+                winner = method(rankings, tiebreaker='random')
+                utility_sums[name][n_cands] += utilities.sum(axis=0)[winner]
 
-        # Find the social utility winner and accumulate utilities
-        UW = utility_winner(utilities)
-        utility_sums['UW'][n_cands] += utilities.sum(axis=0)[UW]
+    return utility_sums
 
-        for name, method in rated_methods.items():
-            winner = method(utilities, tiebreaker='random')
-            utility_sums[name][n_cands] += utilities.sum(axis=0)[winner]
 
-        rankings = honest_rankings(utilities)
-        for name, method in ranked_methods.items():
-            winner = method(rankings, tiebreaker='random')
-            utility_sums[name][n_cands] += utilities.sum(axis=0)[winner]
 
-elapsed_time = time.monotonic() - start_time
-print('Elapsed:', time.strftime("%H:%M:%S", time.gmtime(elapsed_time)), '\n')
+jobs = [delayed(simulate_batch)()] * n_batches
+print(f'{len(jobs)} tasks total:')
+results = Parallel(n_jobs=-3, verbose=5, backend='loky')(jobs)
+
+# Aggregate results
+utility_sums = {key: Counter() for key in (ranked_methods.keys() |
+                                           rated_methods.keys() | {'UW'})}
+for result in results:
+    for method, counter in result.items():
+        utility_sums[method].update(counter)
+
 
 plt.figure(f'Effectiveness, {n_voters} voters, {n_elections} elections')
 plt.title('The Effectiveness of Several Voting Systems')
@@ -91,6 +110,20 @@
 print(tabulate(table, 'keys', showindex=n_cands_list,
                tablefmt="pipe", floatfmt='.2f'))
 
+# Reference values from Weber's published table, ordered by n_cands_list.
+# This script reproduces them within tolerance, so test_examples.py checks the
+# computed table against these.  The paper only gives the m -> infinity limit
+# for the last row, so the 255-candidate value is taken from the "Typical
+# result" table in the docstring instead.
+reference_table = {
+    'Standard': (81.65, 75.00, 69.28, 64.55, 60.61, 49.79, 12.78),
+    'Vote-for-half': (81.65, 75.00, 80.00, 79.06, 81.32, 82.99, 86.37),
+    'Borda': (81.65, 86.60, 89.44, 91.29, 92.58, 95.35, 99.80),
+}
+
+# Absolute tolerance (percentage points) for test_examples.py
+tolerance = 5.0
+
 plt.plot([], [], 'k:', lw=0.8, label='Weber')  # Dummy plot for label
 plt.legend()
 plt.grid(True, color='0.7', linestyle='-', which='major', axis='both')
diff --git a/examples/weber_1977_table_4.py b/examples/weber_1977_table_4.py
index c6bcf2d..3937f5c 100644
--- a/examples/weber_1977_table_4.py
+++ b/examples/weber_1977_table_4.py
@@ -25,10 +25,10 @@
 | 30 |    16.0177 | 16.1577 |    16.1669 |
 
 """
-import time
 from collections import Counter
 
 import numpy as np
+from joblib import Parallel, delayed
 from tabulate import tabulate
 
 from elsim.elections import random_utilities
@@ -39,30 +39,49 @@
 n_voters_list = (2, 3, 4, 5, 10, 15, 20, 25, 30)
 n_cands = 3
 
+# Simulate more than just one election per worker to improve efficiency
+batch_size = 100
+n_batches = n_elections // batch_size
+assert n_batches * batch_size == n_elections
+
 ranked_methods = {'Standard': fptp, 'Borda': borda}
 rated_methods = {'Approval': lambda utilities, tiebreaker:
                  approval(approval_optimal(utilities), tiebreaker)}
 
-utility_sums = {key: Counter() for key in (ranked_methods.keys() |
-                                           rated_methods.keys())}
 
-start_time = time.monotonic()
+def simulate_batch():
+    """Run one batch of elections and return the partial tallies."""
+    utility_sums = {key: Counter() for key in (ranked_methods.keys() |
+                                               rated_methods.keys())}
+
+    for _iteration in range(batch_size):
+        for n_voters in n_voters_list:
+            utilities = random_utilities(n_voters, n_cands)
+
+            for name, method in rated_methods.items():
+                winner = method(utilities, tiebreaker='random')
+                utility_sums[name][n_voters] += utilities.sum(axis=0)[winner]
+
+            rankings = honest_rankings(utilities)
+            for name, method in ranked_methods.items():
+                winner = method(rankings, tiebreaker='random')
+                utility_sums[name][n_voters] += utilities.sum(axis=0)[winner]
 
-for _iteration in range(n_elections):
-    for n_voters in n_voters_list:
-        utilities = random_utilities(n_voters, n_cands)
+    return utility_sums
 
-        for name, method in rated_methods.items():
-            winner = method(utilities, tiebreaker='random')
-            utility_sums[name][n_voters] += utilities.sum(axis=0)[winner]
 
-        rankings = honest_rankings(utilities)
-        for name, method in ranked_methods.items():
-            winner = method(rankings, tiebreaker='random')
-            utility_sums[name][n_voters] += utilities.sum(axis=0)[winner]
 
-elapsed_time = time.monotonic() - start_time
-print('Elapsed:', time.strftime("%H:%M:%S", time.gmtime(elapsed_time)), '\n')
+jobs = [delayed(simulate_batch)()] * n_batches
+print(f'{len(jobs)} tasks total:')
+results = Parallel(n_jobs=-3, verbose=5, backend='loky')(jobs)
+
+# Aggregate results
+utility_sums = {key: Counter() for key in (ranked_methods.keys() |
+                                           rated_methods.keys())}
+for result in results:
+    for method, counter in result.items():
+        utility_sums[method].update(counter)
+
 
 table = {}
 
@@ -73,3 +92,18 @@
 
 print(tabulate(table, 'keys', showindex=n_voters_list,
                tablefmt="pipe", floatfmt='.4f'))
+
+# Reference values from Weber's published Table 4, ordered by n_voters_list.
+# This script reproduces them within tolerance, so test_examples.py checks the
+# computed table against these.
+reference_table = {
+    'Standard': (1.2500, 1.8333, 2.3889, 2.9167, 5.5975, 8.2245,
+                 10.8328, 13.4328, 16.0190),
+    'Borda': (1.2917, 1.8750, 2.4236, 2.9765, 5.6706, 8.3206,
+              10.9472, 13.5588, 16.1597),
+    'Approval': (1.2917, 1.8646, 2.4213, 2.9726, 5.6719, 8.3245,
+                 10.9531, 13.5662, 16.1684),
+}
+
+# Absolute tolerance (utility units) for test_examples.py
+tolerance = 0.3
diff --git a/tests/test_examples.py b/tests/test_examples.py
index 20e6104..55f8240 100644
--- a/tests/test_examples.py
+++ b/tests/test_examples.py
@@ -5,8 +5,12 @@
 ``table`` is read directly, rather than parsing printed output.  These tests
 are slow, so they are marked "slow" and skipped by default (run with
 ``pytest -m slow``), and they require the optional ``examples`` dependencies
-(joblib and matplotlib).  Reference values and tolerances are described at
-``REFERENCE_VALUES``.
+(joblib and matplotlib).
+
+Each example script that has expected output defines ``reference_table`` (the
+values to check its computed ``table`` against) and ``tolerance`` (the
+absolute tolerance for the comparison), so the script doubles as a test.  See
+issue #91.
 """
 import os
 import pathlib
@@ -38,94 +42,13 @@
     'weber_1977_table_4.py',
 ]
 
-# Reference results for the example scripts.
-#
-# For merrill_1984_table_1, merrill_1984_table_3, weber_1977_effectiveness
-# and weber_1977_table_4 the scripts reproduce the published tables within
-# ~2 pp, so the reference values are taken from the papers.
-#
-# merrill_1984_table_2 and merrill_1984_table_4 do not match the published
-# tables (unresolved discrepancies of up to ~5 pp and ~9 pp), so those two are
-# checked against the "Typical result" tables in their docstrings instead.
-# weber_1977_effectiveness_table's last row (255 candidates) is taken from its
-# docstring too, because the paper only gives the m -> infinity limit there.
-#
-# Values are keyed by method and appear in the same order as the script's
-# columns (n_cands, n_voters, or condition).
-REFERENCE_VALUES = {
-    'merrill_1984_table_1_fig_1.py': {
-        'Plurality': (100.0, 79.1, 69.4, 62.1, 52.0, 42.6),
-        'Runoff': (100.0, 96.2, 90.1, 83.6, 73.5, 61.3),
-        'Hare': (100.0, 96.2, 92.7, 89.1, 84.8, 77.9),
-        'Approval': (100.0, 76.0, 69.8, 67.1, 63.7, 61.3),
-        'Borda': (100.0, 90.8, 87.3, 86.2, 85.3, 84.3),
-        'Coombs': (100.0, 96.3, 93.4, 90.2, 86.1, 81.1),
-        'Black': (100.0, 100.0, 100.0, 100.0, 100.0, 100.0),
-        'SU max': (100.0, 84.4, 80.2, 77.9, 77.2, 77.8),
-        'CW': (100.0, 91.6, 83.4, 75.8, 64.3, 52.5),
-    },
-    'merrill_1984_table_2.py': {
-        'Plurality': (57.5, 65.8, 62.2, 78.4, 21.7, 24.4, 27.2, 41.3),
-        'Runoff': (80.1, 87.3, 81.6, 93.6, 35.4, 42.2, 41.5, 61.5),
-        'Hare': (79.2, 86.7, 84.0, 95.4, 35.9, 46.8, 41.0, 69.9),
-        'Approval': (73.8, 77.8, 76.9, 85.4, 71.5, 76.4, 73.8, 82.7),
-        'Borda': (87.1, 89.3, 88.2, 92.3, 83.7, 86.3, 85.2, 89.4),
-        'Coombs': (97.8, 97.3, 97.9, 98.2, 93.5, 92.3, 93.8, 94.5),
-        'Black': (100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0),
-        'SU max': (82.9, 85.8, 85.3, 90.8, 78.1, 81.5, 80.8, 87.1),
-        'CW': (99.7, 99.7, 99.7, 99.6, 98.9, 98.6, 98.7, 98.5),
-    },
-    'merrill_1984_table_3_fig_3.py': {
-        'Plurality': (100.0, 83.0, 75.0, 69.2, 62.8, 53.3),
-        'Runoff': (100.0, 89.5, 83.8, 80.5, 75.6, 67.6),
-        'Hare': (100.0, 89.5, 84.7, 82.4, 80.5, 74.9),
-        'Approval': (100.0, 95.4, 91.1, 89.1, 87.8, 87.0),
-        'Borda': (100.0, 94.8, 94.1, 94.4, 95.4, 95.9),
-        'Coombs': (100.0, 89.7, 86.7, 85.1, 83.1, 82.4),
-        'Black': (100.0, 93.1, 91.9, 92.0, 93.1, 94.3),
-    },
-    'merrill_1984_table_4.py': {
-        'Plurality': (72.1, 79.1, 80.4, 92.4, 4.0, 6.3, 25.2, 52.9),
-        'Runoff': (90.5, 94.2, 92.0, 97.5, 36.6, 43.6, 53.3, 75.3),
-        'Hare': (91.7, 94.7, 94.3, 98.4, 46.4, 57.7, 58.7, 83.6),
-        'Approval': (96.2, 97.0, 96.8, 98.5, 95.6, 96.8, 95.8, 98.0),
-        'Borda': (97.8, 98.6, 98.3, 99.4, 96.6, 97.7, 97.4, 99.0),
-        'Coombs': (97.0, 97.5, 97.7, 98.7, 94.0, 94.3, 95.0, 96.7),
-        'Black': (97.3, 97.8, 98.0, 99.0, 95.5, 96.1, 96.5, 98.0),
-    },
-    'weber_1977_effectiveness_table.py': {
-        'Standard': (81.65, 75.00, 69.28, 64.55, 60.61, 49.79, 12.78),
-        'Vote-for-half': (81.65, 75.00, 80.00, 79.06, 81.32, 82.99, 86.37),
-        'Borda': (81.65, 86.60, 89.44, 91.29, 92.58, 95.35, 99.80),
-    },
-    'weber_1977_table_4.py': {
-        'Standard': (1.2500, 1.8333, 2.3889, 2.9167, 5.5975, 8.2245,
-                     10.8328, 13.4328, 16.0190),
-        'Borda': (1.2917, 1.8750, 2.4236, 2.9765, 5.6706, 8.3206,
-                  10.9472, 13.5588, 16.1597),
-        'Approval': (1.2917, 1.8646, 2.4213, 2.9726, 5.6719, 8.3245,
-                     10.9531, 13.5662, 16.1684),
-    },
-}
-
-# Absolute tolerance per script: percentage points for the Merrill and Weber
-# effectiveness tables, utility units for weber_1977_table_4.
-TOLERANCES = {
-    'merrill_1984_table_1_fig_1.py': 4.0,
-    'merrill_1984_table_2.py': 3.5,
-    'merrill_1984_table_3_fig_3.py': 3.5,
-    'merrill_1984_table_4.py': 3.5,
-    'weber_1977_effectiveness_table.py': 5.0,
-    'weber_1977_table_4.py': 0.3,
-}
-
-
 def _run(name, tmp_path):
-    """Run an example script in full (subprocess); return its ``table``."""
+    """Run an example script in full (subprocess); return its globals."""
     script = (EXAMPLES / name).read_text()
-    out = tmp_path / 'table.pkl'
+    out = tmp_path / 'result.pkl'
     script += (f'\nimport pickle\n'
-               f'pickle.dump(table, open({str(out)!r}, "wb"))\n')
+               f'pickle.dump((table, reference_table, tolerance), '
+               f'open({str(out)!r}, "wb"))\n')
     variant = tmp_path / name
     variant.write_text(script)
     env = {**os.environ, 'MPLBACKEND': 'Agg', 'PYTHONPATH': str(EXAMPLES)}
@@ -135,21 +58,65 @@ def _run(name, tmp_path):
         return pickle.load(f)
 
 
-def _table_rows(table):
-    """Return an example script's ``table`` as {label: np.array}."""
-    if isinstance(table, dict):
-        return {k: np.asarray(v, dtype=float) for k, v in table.items()}
-    return {row[0]: np.asarray(row[1:], dtype=float) for row in table}
+def _method_values(rows):
+    """Return an example script's ``table``/``reference_table`` rows as
+    {method: np.array}.
+
+    ``rows`` is either a list of [method, *values] rows (Merrill-style
+    tables) or a dict mapping method to values (Weber-style tables).
+    """
+    if isinstance(rows, dict):
+        return {k: np.asarray(v, dtype=float) for k, v in rows.items()}
+    return {row[0]: np.asarray(row[1:], dtype=float) for row in rows}
+
+
+def _reference_values(reference):
+    """Return an example script's ``reference_table`` as {method: np.array}.
+
+    Values may be plain sequences in column order, or dicts keyed by column
+    label (sorted by key into column order, e.g. ``merrill_table_1``).
+    """
+    out = {}
+    for method, values in reference.items():
+        if isinstance(values, dict):
+            out[method] = np.asarray(
+                [v for _, v in sorted(values.items())], dtype=float)
+        else:
+            out[method] = np.asarray(values, dtype=float)
+    return out
+
+
+def _is_nested(table):
+    """True if ``table`` is a figure script's dict of {fig: [rows]}."""
+    return isinstance(table, dict) and all(
+        isinstance(v, list) for v in table.values())
+
+
+def _assert_close(name, got, expected, tolerance):
+    """Check one computed row against its reference and report clear errors."""
+    assert got, f'{name}: produced an empty table'
+    for method, expected_row in expected.items():
+        assert method in got, (
+            f'{name}: computed table is missing row {method!r}')
+        assert len(got[method]) == len(expected_row), (
+            f'{name}: row {method!r} has {len(got[method])} values, '
+            f'expected {len(expected_row)}')
+        np.testing.assert_allclose(got[method], expected_row,
+                                   atol=tolerance)
 
 
 @pytest.mark.slow
 @pytest.mark.parametrize('name', ALL_SCRIPTS)
 def test_example(name, tmp_path):
-    table = _table_rows(_run(name, tmp_path))
-    assert table, f'{name}: produced an empty table'
-    for method, expected in REFERENCE_VALUES.get(name, {}).items():
-        assert len(table[method]) == len(expected), (
-            f'{name}: row {method!r} has {len(table[method])} values, '
-            f'expected {len(expected)}')
-        np.testing.assert_allclose(table[method], expected,
-                                   atol=TOLERANCES[name])
+    """Run each example script and check its ``table`` against the
+    ``reference_table``/``tolerance`` defined in the script itself."""
+    table, reference_table, tolerance = _run(name, tmp_path)
+    if _is_nested(table):
+        for fig, rows in table.items():
+            got = _method_values(rows)
+            expected = _reference_values(reference_table[fig])
+            _assert_close(f'{name} ({fig})', got, expected, tolerance)
+    else:
+        got = _method_values(table)
+        expected = _reference_values(reference_table)
+        _assert_close(name, got, expected, tolerance)