#!/usr/bin/python
#
# Copyright (c) 2013 Peter Kasson
# Based on theory in Kasson, Ensign, and Pande, JACS 131(32), 11338-11340.
"""Tool to estimate rate based on single-event observations."""
import math
from scipy import inf
from scipy import integrate
def p_rate(k, go_sequence, stay_sequence):
"""Calculate probability of a rate k given single-event observations.
Args:
k: rate constant to test
go_sequence: list of transition times
stay_sequence: list of observation times with no transition.
Returns:
prob: probability of rate k
"""
theta = sum(go_sequence) + sum(stay_sequence)
n = len(go_sequence)
prob = theta ** (n + 1) /math.factorial(n) * k ** n * math.exp(-k * theta)
return prob
def p_rate_above(k, go_sequence, stay_sequence):
"""Calculate probability of a rate > k given single-event observations.
Args:
k: rate constant to test
go_sequence: list of transition times
stay_sequence: list of observation times with no transition.
Returns:
prob: probability of rate
"""
return integrate.quad(p_rate, k, inf, args = (go_sequence, stay_sequence))[0]
def p_rate_below(k, go_sequence, stay_sequence):
"""Calculate probability of a rate < k given single-event observations.
Args:
k: rate constant to test
go_sequence: list of transition times
stay_sequence: list of observation times with no transition.
Returns:
prob: probability of rate
"""
return integrate.quad(p_rate, 0, k, args = (go_sequence, stay_sequence))[0]