diff --git a/examples/README.md b/examples/README.md
index 9e38b07..2c3c815 100644
--- a/examples/README.md
+++ b/examples/README.md
@@ -2,6 +2,15 @@
 
 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.
+
 ## 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..8f247ab 100644
--- a/examples/merrill_1984_fig_2c_2d.py
+++ b/examples/merrill_1984_fig_2c_2d.py
@@ -87,6 +87,8 @@
     'Plurality': {2: 100.0, 3: 51.3,  4: 36.2,  5: 21.0,  7: 7.8},
 }
 
+table = {}
+
 for fig, disp, ymin, orig in (('2.c', 1.0, 50, merrill_fig_2c),
                               ('2.d', 0.5, 0, merrill_fig_2d)):
 
@@ -134,7 +136,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 +145,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 +158,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..083bddf 100644
--- a/examples/merrill_1984_fig_2c_2d_updated.py
+++ b/examples/merrill_1984_fig_2c_2d_updated.py
@@ -76,6 +76,8 @@
                      star(honest_normed_scores(utilities, 5), tiebreaker),
                  }
 
+table = {}
+
 for fig, disp, ymin in (('2.c', 1.0, 50),
                         ('2.d', 0.5, 0)):
 
@@ -113,7 +115,7 @@
     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 +124,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 +136,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..9cac6f0 100644
--- a/examples/merrill_1984_fig_4a_4b.py
+++ b/examples/merrill_1984_fig_4a_4b.py
@@ -87,6 +87,8 @@
     'Plurality': {2: 100.0, 3: 41.1, 4: 27.0, 5: -1.0, 7: -9},
 }
 
+table = {}
+
 for fig, disp, ymin, orig in (('4.a', 1.0, 55, merrill_fig_4a),
                               ('4.b', 0.5, 0, merrill_fig_4b)):
 
@@ -133,7 +135,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 +145,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 +158,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..4f4e42e 100644
--- a/examples/merrill_1984_fig_4a_4b_updated.py
+++ b/examples/merrill_1984_fig_4a_4b_updated.py
@@ -78,8 +78,10 @@
                           tiebreaker),
                  }
 
+table = {}
+
 for fig, disp, _ymin in (('4.a', 1.0, 55),
-                        ('4.b', 0.5, 0)):
+                         ('4.b', 0.5, 0)):
 
     utility_sums = {key: Counter() for key in (ranked_methods.keys() |
                                                rated_methods.keys() |
@@ -114,7 +116,7 @@
     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 +126,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 +138,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..f35a626 100644
--- a/examples/merrill_1984_table_1_fig_1.py
+++ b/examples/merrill_1984_table_1_fig_1.py
@@ -84,8 +84,10 @@
 elapsed_time = time.monotonic() - start_time
 print('Elapsed:', time.strftime("%H:%M:%S", time.gmtime(elapsed_time)), '\n')
 
-# Plot Merrill's results as dotted lines for comparison
-merrill_table_1 = {
+# 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 +99,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..02eff90 100644
--- a/examples/merrill_1984_table_2.py
+++ b/examples/merrill_1984_table_2.py
@@ -124,3 +124,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..9960447 100644
--- a/examples/merrill_1984_table_3_fig_3.py
+++ b/examples/merrill_1984_table_3_fig_3.py
@@ -79,8 +79,10 @@
 elapsed_time = time.monotonic() - start_time
 print('Elapsed:', time.strftime("%H:%M:%S", time.gmtime(elapsed_time)), '\n')
 
-# Plot Merrill's results as dotted lines for comparison
-merrill_table_1 = {
+# 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 +92,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..e8fcb8a 100644
--- a/examples/merrill_1984_table_4.py
+++ b/examples/merrill_1984_table_4.py
@@ -121,3 +121,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..859d491 100644
--- a/examples/weber_1977_effectiveness_table.py
+++ b/examples/weber_1977_effectiveness_table.py
@@ -91,6 +91,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..1758c8f 100644
--- a/examples/weber_1977_table_4.py
+++ b/examples/weber_1977_table_4.py
@@ -73,3 +73,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 b2dfe4b..ce7bdd9 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,29 +58,138 @@ 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'
+    assert set(got) == set(expected), (
+        f'{name}: computed methods {sorted(got)} do not match reference '
+        f'methods {sorted(expected)}')
+    for method, expected_row in expected.items():
+        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)
+
+
+def _check_nested(name, table, reference_table, tolerance):
+    """Check a figure script's computed ``table`` against its reference."""
+    assert set(table) == set(reference_table), (
+        f'{name}: computed figures {sorted(table)} do not match '
+        f'reference figures {sorted(reference_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)
 
 
 @pytest.mark.slow
 @pytest.mark.parametrize('name', ALL_SCRIPTS)
 def test_example(name, tmp_path):
-    """Run an example script in a subprocess and check its ``table``.
+    """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):
+        _check_nested(name, table, reference_table, tolerance)
+    else:
+        got = _method_values(table)
+        expected = _reference_values(reference_table)
+        _assert_close(name, got, expected, tolerance)
+
+
+def test_assert_close_rejects_extra_method():
+    """A computed method with no reference row must fail the test.
 
-    The script is run to completion and its computed ``table`` is compared
-    against ``REFERENCE_VALUES`` within ``TOLERANCES``.  A non-empty table
-    and equal-length rows are asserted first, so a script that silently
-    truncates its output fails with a clear message rather than an opaque
-    array mismatch.
+    Previously only missing methods were caught; an extra method could be
+    silently ignored.
     """
-    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])
+    got = {'A': np.array([1.0]), 'B': np.array([1.0])}
+    with pytest.raises(AssertionError):
+        _assert_close('x', got, {'A': np.array([1.0])}, 0.1)
+
+
+def test_assert_close_rejects_missing_method():
+    """A reference method absent from the computed table must fail."""
+    got = {'A': np.array([1.0])}
+    with pytest.raises(AssertionError):
+        _assert_close('x', got, {'A': np.array([1.0]),
+                                 'B': np.array([2.0])}, 0.1)
+
+
+def test_assert_close_rejects_wrong_length():
+    """A row with the wrong number of values must fail with a clear message."""
+    got = {'A': np.array([1.0])}
+    with pytest.raises(AssertionError, match='has 1 values'):
+        _assert_close('x', got, {'A': np.array([1.0, 2.0])}, 0.1)
+
+
+def test_assert_close_rejects_different_values():
+    """Values differing beyond ``tolerance`` must fail."""
+    got = {'A': np.array([1.0])}
+    with pytest.raises(AssertionError):
+        _assert_close('x', got, {'A': np.array([2.0])}, 0.1)
+
+
+def test_assert_close_rejects_empty_table():
+    """An empty computed table must fail."""
+    with pytest.raises(AssertionError, match='empty table'):
+        _assert_close('x', {}, {'A': np.array([1.0])}, 0.1)
+
+
+def test_check_nested_rejects_missing_figure():
+    """A reference figure absent from the computed table must fail.
+
+    Without this check a missing figure would be silently ignored.
+    """
+    table = {'2.c': [['A', 1.0]]}
+    reference = {'2.c': {'A': np.array([1.0])},
+                 '2.d': {'A': np.array([1.0])}}
+    with pytest.raises(AssertionError, match='computed figures'):
+        _check_nested('x', table, reference, 0.1)
+
+
+def test_check_nested_rejects_empty_nested_table():
+    """An empty nested table must fail rather than pass vacuously."""
+    with pytest.raises(AssertionError, match='computed figures'):
+        _check_nested('x', {}, {'2.c': {'A': np.array([1.0])}}, 0.1)
+
+
+def test_check_nested_passes_matching_figures():
+    """A nested table matching its reference must pass."""
+    table = {'2.c': [['A', 1.0]]}
+    reference = {'2.c': {'A': np.array([1.0])}}
+    _check_nested('x', table, reference, 0.1)