testing/bell_rate/fit_bell

   1 #!/usr/bin/python
   2 #
   3 # fit data coming in on stdin to bell-theory histogram
   4 # output:
   5 #     <convergence information>
   6 #   a:\t<some number>
   7 #   b:\t<another number>
   8 #   N:\t<number of data points>
   9 #   bin:\t<width of histogram bins>
  10 #
  11 # usage: fit_bell
  12 # for example:
  13 #   $ cat datafile | fit_bell
  14
  15 import sys
  16 sys.path.append('../common/')
  17 import fit
  18 from scipy import exp
  19
  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) :
  27         a = params[0]
  28         b = params[1]
  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"
  32         init = self.data[0]
  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])
  36         Ymax = max(Yarr)
  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
  44             i = 0
  45             s = 0
  46             while s < self.N/2.0 :
  47                 s += Yarr[i]
  48                 i += 1
  49             print "half dead by", self.data[i]
  50             Aguess = init[1]/self.N/self.binwidth
  51             Bguess = 1.0/self.data[i][0]
  52         else :
  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
  58             #   b = a e^(b fmax)
  59             # hmm, trancendental equation for b
  60             print "Xmax = ", Xmax
  61             Bguess = 1.0/Xmax
  62             Aguess = init[1]/self.N/self.binwidth
  63         print "guessing ", (Aguess, Bguess)
  64         return (Aguess, Bguess)
  65
  66 if __name__ == "__main__" :
  67     fit.run(2, Model)