from math import log, log1p, exp, trunc, ceil, floor
from numbers import Real, Integral
import numpy as np
[docs]
class Probability(Real):
"""Probability class.
Similar to a float, but value stored as natural log value internally.
All operations are attempted with log values where possible, and failing
that a float will be returned instead.
Parameters
----------
value : float
Value for probability.
log_value : bool
Set to `True` if :attr:`value` already a log value. Default `False`.
"""
def __init__(self, value, *, log_value=False):
if log_value:
self._log_value = value
else:
try:
self._log_value = self._log(value)
except ValueError:
raise ValueError("value must be greater than 0")
@property
def log_value(self):
return self._log_value
@staticmethod
def _log(other):
if isinstance(other, Probability):
return other.log_value
elif other == 0:
return float("-inf")
else:
return log(other)
def __hash__(self):
value = float(self)
if value == 0 and self.log_value != float("-inf"): # Too close to zero
# Add string so doesn't have same hash as the log value itself.
return hash(('log', self._log_value))
else:
return hash(value)
def __eq__(self, other):
if other < 0:
return False
return self.log_value == self._log(other)
def __le__(self, other):
if other < 0:
return False
return self.log_value <= self._log(other)
def __lt__(self, other):
if other < 0:
return False
return self.log_value < self._log(other)
def __ge__(self, other):
if other < 0:
return True
return self.log_value >= self._log(other)
def __gt__(self, other):
if other < 0:
return True
return self.log_value > self._log(other)
def __add__(self, other):
if other < 0:
return self - -other
log_other = self._log(other)
if self.log_value > log_other:
log_l, log_s = self.log_value, log_other
elif self.log_value < log_other:
log_l, log_s = log_other, self.log_value
else: # Must be equal, so just double value
return self * 2
if log_s == float("-inf"): # Just return largest value
return Probability(log_l, log_value=True)
return Probability(log_l + log1p(exp(log_s - log_l)),
log_value=True)
def __radd__(self, other):
return self + other
def __sub__(self, other):
if other < 0:
return self + -other
log_other = self._log(other)
if self.log_value > log_other:
log_l, log_s = self.log_value, log_other
elif self.log_value < log_other: # Result will be negative
return float(self) - other
else: # Must be equal, so return 0
return Probability(float("-inf"), log_value=True)
if log_s == float("-inf"): # Just return largest value
return Probability(log_l, log_value=True)
exp_diff = exp(log_s - log_l)
if exp_diff == 1: # Diff too small, so result is effectively zero
return Probability(float("-inf"), log_value=True)
return Probability(log_l + log1p(-exp_diff),
log_value=True)
def __rsub__(self, other):
if other < 0: # Result will be negative
return other + -float(self)
log_other = self._log(other)
if log_other > self.log_value:
log_l, log_s = log_other, self.log_value
elif log_other < self.log_value: # Result will be negative
return other + -float(self)
else: # Must be equal, so return 0
return Probability(float("-inf"), log_value=True)
if log_s == float("-inf"): # Just return largest value
return Probability(log_l, log_value=True)
exp_diff = exp(log_s - log_l)
if exp_diff == 1: # Diff too small, so result is effectively zero
return Probability(float("-inf"), log_value=True)
return Probability(log_l + log1p(-exp_diff),
log_value=True)
def __imul__(self, other):
try:
return Probability(self.log_value + self._log(other),
log_value=True)
except ValueError:
return float(self) * other
def __mul__(self, other):
try:
if isinstance(other, Probability):
return Probability(self.log_value + self._log(other),
log_value=True)
else:
return exp(self.log_value + self._log(other))
except ValueError:
return float(self) * other
def __rmul__(self, other):
return self * other
def __truediv__(self, other):
try:
return Probability(self.log_value - self._log(other),
log_value=True)
except ValueError:
return float(self) / other
def __rtruediv__(self, other):
try:
return Probability(self._log(other) - self.log_value,
log_value=True)
except ValueError:
return other / float(self)
def __floordiv__(self, other):
return floor(self / other)
def __rfloordiv__(self, other):
return floor(other / self)
def __mod__(self, other):
try:
return Probability(float(self) % other)
except ValueError:
return float(self) % other
def __rmod__(self, other):
return Probability(other % float(self))
def __pow__(self, exponent):
if isinstance(exponent, Probability):
exponent = float(exponent)
return Probability(exponent * self.log_value, log_value=True)
def __rpow__(self, base):
return Probability(base ** float(self))
def __neg__(self):
return -float(self)
def __pos__(self):
return Probability(self)
def __abs__(self):
return Probability(self)
def __float__(self):
return exp(self.log_value)
def __round__(self, ndigits=None):
value = round(float(self), ndigits)
if isinstance(value, Integral):
return value
else:
return Probability(value)
def __trunc__(self):
return trunc(float(self))
def __floor__(self):
return floor(float(self))
def __ceil__(self):
return ceil(float(self))
def __repr__(self):
value = float(self)
if value == 0 and self.log_value != float("-inf"): # Too close to zero
return "Probability({!r}, log_value=True)".format(self.log_value)
else:
return "Probability({!r})".format(value)
def __str__(self):
value = float(self)
if value == 0 and self.log_value != float("-inf"): # Too close to zero
return "exp({})".format(self.log_value)
else:
return str(value)
[docs]
def sqrt(self):
"""Square root which can be called by NumPy"""
return self ** 0.5
[docs]
def log(self):
"""Log which can be called by NumPy"""
return self.log_value
[docs]
@classmethod
def sum(cls, values):
"""Carry out LogSumExp"""
log_values = np.array([cls._log(value) for value in values])
if len(log_values) == 0:
return Probability(0)
max_log_value = max(log_values)
# Check if all values are zero
if max_log_value == float('-inf'):
return Probability(0)
value_sum = np.sum(np.exp(log_values - max_log_value))
return Probability(cls._log(value_sum) + max_log_value, log_value=True)
[docs]
@classmethod
def from_log(cls, value):
"""Create Probability type from a log value; same as Probability(value, log_value=True)
Parameters
----------
value : float
Value for probability.
Returns
-------
Probability
"""
return cls(value, log_value=True)
Probability.from_log_ufunc = np.frompyfunc(Probability.from_log, 1, 1)