In addition to the existing answers, this can actually also be done very easily using matplotlib.pyplot.hexbin.

Demonstration with the same setup as @Andrea:

import matplotlib.pyplot as plt
import numpy as np

x = np.random.uniform(0, 10, 1000)
y = np.random.uniform(10, 20, 1000)
z = np.exp(-(x-3)**2/5 - (y-18)**2/5) + np.random.random(1000)

fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(12, 4))
ax0.scatter(x, y, c=z)
ax1.hexbin(x, y, C=z, gridsize=10)

That last line is all it takes. If keyword argument C is None (the default), then hexbin indeed behaves as a histogram that counts the number of points within each hexagon. However, if C is provided, it does exactly what OP is asking, combining the values within each bin based on reduce_C_function (which by default is the mean).