Credible interval of n dimensional, unimodal distributions in Python

credible_interval_nd.py

import itertools as it

def credible_interval(lims, ys, mass=0.95):
    dimensions = ys.shape
    coordinates = it.product(*(range(x) for x in dimensions))
    indexed = zip(ys.ravel(), coordinates)
    indexed.sort(reverse=True)
    
    acc = 0
    points = []
    for (m, c) in indexed:
        if acc + m > mass:
            break
        acc += m
        points.append(c)

        
    mins = []
    maxs = []
    for (pos, lim, dim) in zip(range(len(dimensions)), lims, dimensions):
        mini = min(points, key=lambda(t): t[pos])[pos]
        maxi = max(points, key=lambda(t): t[pos])[pos]
        span = np.linspace(lim[0], lim[1], dim)
        mins.append(span[mini])
        maxs.append(span[maxi])
    
    return (mins, maxs)