Compare theory vs. sim scaled by bin width (not theory, which may be zero).

[sawsim.git] / pysawsim / test / bell_rate.py

diff --git a/pysawsim/test/bell_rate.py b/pysawsim/test/bell_rate.py
index b7eb3f63c32618644a963bd05713a1a9b0b61d18..b002a9a2d7d9c49be12012dd71467151905e50f7 100644 (file)
--- a/pysawsim/test/bell_rate.py
+++ b/pysawsim/test/bell_rate.py
@@ -96,11 +96,20 @@ def probability_distribution(x, params):
     """Gumbel (minimum) probability distribution.
 
     .. math:: 1/\\rho exp((x-\\alpha)/\\rho - exp((x-\\alpha)/\\rho))
+
+    Integral
+
+    .. math:: -exp(-exp((x-\\alpha)/\\rho))
+
+    So integrated over the range x = [0,\\infty]
+
+    .. math:: -exp(-\\infty) - (-exp(-exp(-\\alpha/\\rho)))
+      = exp(-exp(-\\alpha/\\rho)))
     """
     p = params  # convenient alias
     p[1] = abs(p[1])  # cannot normalize negative rho.
     xs = (x - p[0]) / p[1]
-    return (1.0/p[1]) * exp(xs - exp(xs))
+    return (exp(exp(-p[0]/p[1]))/p[1]) * exp(xs - exp(xs))
 
 
 class GumbelModelFitter (HistogramModelFitter):
@@ -146,7 +155,7 @@ def bell_rate(sawsim_runner, num_domains=1, unfolding_rate=1,
     w = 0.1 * rho # calculate bin width (in force)
     force_mean = alpha - gamma_e * rho
     theory = Histogram()
-    theory.bin_edges = arange(start=0, stop=3*force_mean, step=w)
+    theory.bin_edges = arange(start=0, stop=max(force_mean,0)+3*rho, step=w)
     theory.bin_centers = theory.bin_edges[:-1] + w/2
     theory.counts = w*num_domains*N*probability_distribution(
         theory.bin_centers, [alpha, rho])
@@ -164,13 +173,14 @@ def bell_rate(sawsim_runner, num_domains=1, unfolding_rate=1,
 
     sim = sawsim_histogram(sawsim_runner, param_string, N=N,
                            bin_edges=theory.bin_edges)
-    
+
     e = GumbelModelFitter(sim)
     params = e.fit()
     sim_alpha = params[0]
     sim_rho = abs(params[1])
     for s,t,n in [(sim_alpha, alpha, 'alpha'), (sim_rho, rho, 'rho')]:
-        assert (s - t)/t < 0.1, 'simulation %s = %g != %g = %s' % (n,s,t,n)
+        assert (s - t)/w < 2, (
+            'simulation %s = %g != %g = %s (bin width = %g)' % (n,s,t,n,w))
     return sim.residual(theory)