projects
/
sawsim.git
/ commitdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
| commitdiff |
tree
raw
|
patch
|
inline
| side by side (parent:
bbda76f
)
Fix Gumbel normalization in pysawsim.test.bell_rate.probability_distribution.
author
W. Trevor King
<wking@drexel.edu>
Sat, 6 Nov 2010 14:15:48 +0000
(07:15 -0700)
committer
W. Trevor King
<wking@drexel.edu>
Sat, 6 Nov 2010 14:15:48 +0000
(07:15 -0700)
pysawsim/test/bell_rate.py
patch
|
blob
|
history
diff --git
a/pysawsim/test/bell_rate.py
b/pysawsim/test/bell_rate.py
index 83ca86d9a855f6405024ceafaa798fe95072f84d..a49da052de007a3ddffb14a96763154015af5c52 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))
"""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]
"""
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):
class GumbelModelFitter (HistogramModelFitter):
@@
-164,13
+173,13
@@
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)
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')]:
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)/t < 0.1, 'simulation %s = %g != %g = %s' % (n,s,t,n)
)
return sim.residual(theory)
return sim.residual(theory)