code: def hpd(x, alpha):
      
  
      # Make a copy of trace
      x = x.copy()
  
      # For multivariate node
      if x.ndim > 1:
  
          # Transpose first, then sort
          tx = tr(x, list(range(x.ndim)[1:]) + [0])
          dims = shape(tx)
  
          # Container list for intervals
          intervals = np.resize(0.0, (2,) + dims[:-1])
  
          for index in make_indices(dims[:-1]):
  
              try:
                  index = tuple(index)
              except TypeError:
                  pass
  
              # Sort trace
              sx = sort(tx[index])
  
              # Append to list
              intervals[0][index], intervals[1][index] = calc_min_interval(sx, alpha)
          
          # Transpose back before returning
          return array(intervals)
  
      else:
          # Sort univariate node
          sx = sort(x)
  
          return array(calc_min_interval(sx, alpha))
docs: Calculate highest posterior density (HPD) of array for given alpha. The HPD is the
      minimum width Bayesian credible interval (BCI).
  
      :Arguments:
        x : Numpy array
            An array containing MCMC samples
        alpha : float
            Desired probability of type I error