3 # fit data coming in on stdin to bell-theory histogram
5 # <convergence information>
8 # N:\t<number of data points>
9 # bin:\t<width of histogram bins>
13 # $ cat datafile | fit_bell
16 sys.path.append('../common/')
20 class Model (fit.Model) :
21 def printModel(self) :
22 print "bin*N*a*exp(b*x + a/b(1-exp(b*x)))"
23 def printParams(self, xopt) :
24 print "a:\t%g" % xopt[0]
25 print "b:\t%g" % xopt[1]
26 def model(self, x, params) :
29 return self.binwidth*self.N*a*exp(b*x + a/b*(1.-exp(b*x)))
30 def guessInitialParams(self) :
31 print "guessing initial params"
33 # find the coords (Xmax,Ymax) of the peak of the histogram
34 Yarr = list(zip(*self.data)[1])
35 Xarr = list(zip(*self.data)[0])
37 Xmax = Xarr[Yarr.index(Ymax)]
38 # two guessing methods for B:
39 if Ymax < 1.3*init[1] :
40 # in the first, we have a high unforced rate, and
41 # the histogram is more-or-less exponential decay.
42 print "guessing with method 1"
43 Aguess = init[1]/self.N/self.binwidth
46 while s < self.N/2.0 :
49 print "half dead by", self.data[i]
50 Aguess = init[1]/self.N/self.binwidth
51 Bguess = 1.0/self.data[i][0]
53 # in the second, we have a low unforced rate, and the histogram has
54 # a well defined 'bump' at some likely unfolding force
55 # h(f) = binwidth N a e^{b x + a/b[1-e^(b f)]}
56 # dh/df = h [b - a e^(b f)]
57 # so dh/df = 0 at f=fmax when
59 # hmm, trancendental equation for b
62 Aguess = init[1]/self.N/self.binwidth
63 print "guessing ", (Aguess, Bguess)
64 return (Aguess, Bguess)
66 if __name__ == "__main__" :