Mean values depending on binning with respect to second variable

Question

I am working with python / numpy. As input data I have a large number of value pairs (x,y). I basically want to plot <y>(x), i.e., the mean value of y for a certain data bin x. At the moment I use a plain for loop to achieve this, which is terribly slow.

# create example data
x = numpy.random.rand(1000)
y = numpy.random.rand(1000)
# set resolution
xbins = 100
# find x bins
H, xedges, yedges = numpy.histogram2d(x, y, bins=(xbins,xbins) )
# calculate mean and std of y for each x bin
mean = numpy.zeros(xbins)
std = numpy.zeros(xbins)
for i in numpy.arange(xbins):
    mean[i] = numpy.mean(y[ numpy.logical_and( x>=xedges[i], x<xedges[i+1] ) ])
    std[i]  = numpy.std (y[ numpy.logical_and( x>=xedges[i], x<xedges[i+1] ) ])

Is it possible to have a kind of vectorized writing for it?

Answer 1

You are complicating things unnecessarily. All you need to know is, for every bin in x, what are n, sy and sy2, the number of y values in that x bin, the sum of those y values, and the sum of their squares. You can get those as:

>>> n, _ = np.histogram(x, bins=xbins)
>>> sy, _ = np.histogram(x, bins=xbins, weights=y)
>>> sy2, _ = np.histogram(x, bins=xbins, weights=y*y)

From those:

>>> mean = sy / n
>>> std = np.sqrt(sy2/n - mean*mean)

Answer 2

If you can use pandas:

import pandas as pd
xedges = np.linspace(x.min(), x.max(), xbins+1)
xedges[0] -= 0.00001
xedges[-1] += 0.000001
c = pd.cut(x, xedges)
g = pd.groupby(pd.Series(y), c.labels)
mean2 = g.mean()
std2 = g.std(0)