Changed formula

A slightly different expression is used as it gives slightly better numerical accuracy at the extremes (q=0 and q=1).

Author: vetlewi

ompy/dist/fermi_dirac.py

@@ -53,8 +53,7 @@ def _ppf(self, q, lam, mu):
         """
         Calculate the CDF for the Fermi-Dirac distribution.
         """
-        N = np.log(1 + np.exp(lam*mu))
-        return mu - np.log(np.exp(N*(1-q)) - 1)/lam
+        return mu - np.log((1 + np.exp(lam*mu))**(1-q) - 1)/lam
 
     def _random(self, lam, mu, size=None):
         """
@@ -101,6 +100,10 @@ def logp(self, value):
         TensorVariable
         """
 
+        # The formula is
+        # p(x) = lam/ln(1 + exp(lam*mu)) * 1/(exp(lam*(x-mu)) + 1)
+        # ln(p(x)) = ln(lam/ln(1 + exp(lam*mu))) - ln(exp(lam*(x-mu)) + 1)
+
         lam = self.lam
         mu = self.mu
 
@@ -125,6 +128,10 @@ def logcdf(self, value):
         TensorVariable
         """
 
+        # The formula is CDF
+        # P(x) = 1 - ln(1 + exp(-lam*(x-mu)))/ln(1 + exp(lam*mu))
+        # ln(P(x)) = ln(1 - ln(1 + exp(-lam*(x-mu)))/ln(1 + exp(lam*mu)))
+
         lam = self.lam
         mu = self.mu
 

ompy/dist/fermi_dirac_pymc4.py

@@ -37,8 +37,7 @@ def rng_fn(cls,
                mu: np.ndarray, size: Tuple[int, ...]
                ) -> np.ndarray:
         q = rng.uniform(size=size)
-        N = lam/np.log(1 + np.exp(lam*mu))
-        return mu - np.log(1 - np.exp(lam*(1 - q)/lam))/lam
+        return mu - np.log((1 + np.exp(lam*mu))**(1-q) - 1)/lam
 
 
 fermidirac = FermiDiracRV()
@@ -122,7 +121,6 @@ def logcdf(value, lam, mu):
         TensorVariable
         """
 
-        N = lam/at.log(1 + at.exp(lam*mu))
-        V = (at.exp(lam*(value - mu)) + 1)/(at.exp(lam*mu) + 1)
-        logcdf = at.log(N) + at.log(lam*value - at.log(V))
+        logcdf = at.log(1 - at.log(1 + at.exp(-lam*(value - mu))) /
+                        at.log(1 + at.exp(lam*mu)))
         return bound(logcdf, value >= 0, lam > 0)

tests/test_dist.py

@@ -6,6 +6,7 @@
 
 import matplotlib.pyplot as plt
 
+
 @pytest.mark.parametrize(
     "lam,mu",
     [(10., 1.),
@@ -29,6 +30,7 @@ def test_fermi_dirac(lam, mu):
     # We use a fairly large tolerance since error should be stocastic
     assert_allclose(hist, y, atol=0.1)
 
+
 @pytest.mark.parametrize(
     "lam,mu",
     [(10., 1.),