I went ahead and implemented a solution to this, doing what I thought was exactly what you suggested. And to my surprise, it turns out to execute 6 times slower than my original implementation! i.e. This code:

```
    ecenSM = get_bin_centers(ebinsSM)
    czcenSM = get_bin_centers(czbinsSM)
    map_shape = (len(ecenSM),len(czcenSM))
    map_sum_wts = np.zeros(map_shape, np.float32)
    map_num = np.zeros(map_shape, np.float32)
    indexi = np.digitize(ecenLT,ebinsSM) - 1
    indexj = np.digitize(czcenLT,czbinsSM) - 1

    for ie,i in enumerate(indexi):
        for icz,j in enumerate(indexj):
            map_sum_wts[i,j] += prob_map[ie,icz]
            map_num[i,j] += 1

   return np.divide(map_sum_wts,map_num)
```

is 6 times slower than this:

```
    elist = []                                                                           
    czlist = []                                                                          
    weight_list = []                                                                     
    for ie,egy in enumerate(ecenLT):                                                     
        for icz,cz in enumerate(czcenLT):                                                
            czlist.append(cz)                                                            
            elist.append(egy)                                                            
            weight_list.append(prob_map[ie][icz])                                        

    map_sum_wts = np.histogram2d(elist,czlist,weights=weight_list,                       
                                 bins=[ebinsSM,czbinsSM])[0]                             
    map_num = np.histogram2d(elist,czlist,bins=[ebinsSM,czbinsSM])[0]                    
    return np.divide(map_sum_wts,map_num)
```
Do you have any idea why this would be? Did you have another implementation in mind for the "np.digitize() way"?