In [ ]:

    

import numpy as np


    

In [4]:

    

def fastmultihist(sample, bins):
    N, D = sample.shape
    # Compute the bin number each sample falls into.
    Ncount = {}
    for i in arange(D):
        Ncount[i] = np.digitize(sample[:, i], bins[i])

    for i in arange(D):
        nbin[i] = len(bins[i]) - 1
        
    xy = zeros(N, int)
    for i in arange(0, D-1):
        xy += Ncount[ni[i]] * nbin[ni[i+1:]].prod()
    xy += Ncount[ni[-1]]
    
    hist = zeros(nbin, float).reshape(-1)
    flatcount = bincount(xy, weights)
    a = arange(len(flatcount))
    hist[a] = flatcount
    
    hist = hist.reshape(nbin)


    

In [ ]: