diff --git a/examples/merrill_1984_fig_2c_2d.py b/examples/merrill_1984_fig_2c_2d.py
index 114586c..98bf8cc 100644
--- a/examples/merrill_1984_fig_2c_2d.py
+++ b/examples/merrill_1984_fig_2c_2d.py
@@ -44,6 +44,7 @@
 import time
 from collections import Counter
 import numpy as np
+from scipy.stats import binomtest
 import matplotlib.pyplot as plt
 from tabulate import tabulate
 from elsim.methods import (fptp, runoff, irv, approval, borda, coombs,
@@ -140,7 +141,14 @@
                    'Plurality'):
         x, y = zip(*sorted(condorcet_winner_count[method].items()))
         CE = np.array(y)/y_cw
-        plt.plot(x, CE*100, '-', label=method)
+
+        # Add 95% confidence interval error bars (Clopper-Pearson exact method)
+        ci = np.empty((2, len(y)))
+        for i in range(len(y)):
+            ci[:, i] = binomtest(y[i], y_cw[i]).proportion_ci()
+        yerr = ci - CE
+        yerr[0] = -yerr[0]
+        plt.errorbar(x, CE*100, yerr*100, fmt='-', label=method)
         table.append([method, *CE*100])
 
     print(tabulate(table, ["Method", *x], tablefmt="pipe", floatfmt='.1f'))
diff --git a/examples/merrill_1984_table_1_fig_1.py b/examples/merrill_1984_table_1_fig_1.py
index db983bb..8861ff1 100644
--- a/examples/merrill_1984_table_1_fig_1.py
+++ b/examples/merrill_1984_table_1_fig_1.py
@@ -27,6 +27,7 @@
 import time
 from collections import Counter
 import numpy as np
+from scipy.stats import binomtest
 import matplotlib.pyplot as plt
 from tabulate import tabulate
 from elsim.methods import (fptp, runoff, irv, approval, borda, coombs,
@@ -113,7 +114,14 @@
                'Black'):
     x, y = zip(*sorted(condorcet_winner_count[method].items()))
     CE = np.array(y)/y_cw
-    plt.plot(x, CE*100, '-', label=method)
+
+    # Add 95% confidence interval error bars (Clopper-Pearson exact method)
+    ci = np.empty((2, len(y)))
+    for i in range(len(y)):
+        ci[:, i] = binomtest(y[i], y_cw[i]).proportion_ci()
+    yerr = ci - CE
+    yerr[0] = -yerr[0]
+    plt.errorbar(x, CE*100, yerr*100, fmt='-', label=method)
     table.append([method, *np.array(y)/y_cw*100])
 
 # Likelihood that social utility maximizer is Condorcet Winner
diff --git a/examples/results/Merrill_1984_Figure_1._25_voters,_10000_iterations.png b/examples/results/Merrill_1984_Figure_1._25_voters,_10000_iterations.png
index f0f2aa6..6f3e7b6 100644
Binary files a/examples/results/Merrill_1984_Figure_1._25_voters,_10000_iterations.png and b/examples/results/Merrill_1984_Figure_1._25_voters,_10000_iterations.png differ
diff --git a/examples/results/Merrill_1984_Figure_2.c._201_voters,_10000_iterations.png b/examples/results/Merrill_1984_Figure_2.c._201_voters,_10000_iterations.png
index 8b0c691..39371ba 100644
Binary files a/examples/results/Merrill_1984_Figure_2.c._201_voters,_10000_iterations.png and b/examples/results/Merrill_1984_Figure_2.c._201_voters,_10000_iterations.png differ
diff --git a/examples/results/Merrill_1984_Figure_2.d._201_voters,_10000_iterations.png b/examples/results/Merrill_1984_Figure_2.d._201_voters,_10000_iterations.png
index 994b516..eeb8979 100644
Binary files a/examples/results/Merrill_1984_Figure_2.d._201_voters,_10000_iterations.png and b/examples/results/Merrill_1984_Figure_2.d._201_voters,_10000_iterations.png differ
diff --git a/examples/results/Wikipedia_Condorcet_paradox_likelihood.png b/examples/results/Wikipedia_Condorcet_paradox_likelihood.png
index 96c2b22..95835fd 100644
Binary files a/examples/results/Wikipedia_Condorcet_paradox_likelihood.png and b/examples/results/Wikipedia_Condorcet_paradox_likelihood.png differ
diff --git a/examples/wikipedia_condorcet_paradox_likelihood.py b/examples/wikipedia_condorcet_paradox_likelihood.py
index 3966980..3518858 100644
--- a/examples/wikipedia_condorcet_paradox_likelihood.py
+++ b/examples/wikipedia_condorcet_paradox_likelihood.py
@@ -22,6 +22,7 @@
 """
 from collections import Counter
 import numpy as np
+from scipy.stats import binomtest
 import matplotlib.pyplot as plt
 from tabulate import tabulate
 from joblib import Parallel, delayed
@@ -71,7 +72,14 @@ def simulate_batch(n_voters):
 
 x, y = zip(*sorted(condorcet_paradox_counts.items()))
 CP = np.asarray(y) / n_elections  # Likelihood of paradox
-plt.plot(x, CP*100, '-', label='Simulation')
+
+# Add 95% confidence interval error bars (Clopper-Pearson exact method)
+ci = np.empty((2, len(y)))
+for i in range(len(y)):
+    ci[:, i] = binomtest(y[i], n_elections).proportion_ci()
+yerr = ci - CP
+yerr[0] = -yerr[0]
+plt.errorbar(x, CP*100, yerr*100, fmt='-', label='Simulation')
 
 plt.legend()
 plt.grid(True, color='0.7', linestyle='-', which='major', axis='both')