Redo Fig. 1b

lcs/utilities.py

@@ -1,6 +1,7 @@
 import random
 
 import numpy as np
+from scipy.stats import gaussian_kde
 
 from .contagion import *
 from .inference import *
@@ -67,3 +68,59 @@ def degrees(A):
     if not isinstance(A, np.ndarray):
         A = A.todense()
     return A.sum(axis=0)
+
+
+def hpd_grid(sample, alpha=0.05, roundto=2):
+    """Calculate highest posterior density (HPD) of array for given alpha. 
+    The HPD is the minimum width Bayesian credible interval (BCI). 
+    The function works for multimodal distributions, returning more than one mode
+
+    Parameters
+    ----------
+    
+    sample : Numpy array or python list
+        An array containing MCMC samples
+    alpha : float
+        Desired probability of type I error (defaults to 0.05)
+    roundto: integer
+        Number of digits after the decimal point for the results
+
+    Returns
+    -------
+    hpd: array with the lower 
+
+    References
+    ----------
+    Bayesian Analysis with Python (Second edition)
+    https://github.com/aloctavodia/BAP/blob/master/first_edition/code/Chp1/hpd.py
+    """
+    sample = np.asarray(sample)
+    sample = sample[~np.isnan(sample)]
+    # get upper and lower bounds
+    l = np.min(sample)
+    u = np.max(sample)
+    density = gaussian_kde(sample)
+    x = np.linspace(l, u, 2000)
+    y = density.evaluate(x)
+    #y = density.evaluate(x, l, u) waitting for PR to be accepted
+    xy_zipped = zip(x, y/np.sum(y))
+    xy = sorted(xy_zipped, key=lambda x: x[1], reverse=True)
+    xy_cum_sum = 0
+    hdv = []
+    for val in xy:
+        xy_cum_sum += val[1]
+        hdv.append(val[0])
+        if xy_cum_sum >= (1-alpha):
+            break
+    hdv.sort()
+    diff = (u-l)/20  # differences of 5%
+    hpd = []
+    hpd.append(round(min(hdv), roundto))
+    for i in range(1, len(hdv)):
+        if hdv[i]-hdv[i-1] >= diff:
+            hpd.append(round(hdv[i-1], roundto))
+            hpd.append(round(hdv[i], roundto))
+    hpd.append(round(max(hdv), roundto))
+    ite = iter(hpd)
+    hpd = list(zip(ite, ite))
+    return hpd, x, y