Matplotlib plotting a 3d number density plot using ax.scatter with 2d histograms on axis surfaces?

I have the set of subplots shown below, but I want to apply the 3 ax.hist2d() plots to the 3d projection on the surfaces, is there a way of doing this? My data is a series of coordinates in x, y, z with corresponding kde values calculated from scipy.stats.gaussian_kde() to give the density fluctuations. hist2d has no 'offset' argument for me to add it to the first plot? Do I need to use ax.surface_plot()?

import numpy as np
import matplotlib.pyplot as plt

X = np.random.normal(5, 2, size=(100,))
Y = np.random.normal(5, 2, size=(100,))
Z = np.random.normal(5, 2, size=(100,))
kde  = np.random.normal(1e-6, 1e-7, size=(100,))

fig = plt.figure(figsize=(12,12))

ax = fig.add_subplot(2, 2, 1, projection='3d')
ax.scatter(X, -Z, Y, c=kde, alpha=1, s=1)
ax.invert_xaxis()
ax.set_xlabel('x (cm)')
ax.set_ylabel('z (cm)')
ax.set_zlabel('y (cm)')

ax1 = fig.add_subplot(2, 2, 2)
ax1.hist2d(X, Y, bins=(50,50))
ax1.invert_xaxis()
ax1.set_xlabel('x (cm)')
ax1.set_ylabel('y (cm)')

ax2 = fig.add_subplot(2, 2, 3)
ax2.hist2d(X, -Z, bins=(50,50))
ax2.invert_xaxis()
ax2.set_xlabel('x (cm)')
ax2.set_ylabel('z (cm)')

ax3 = fig.add_subplot(2, 2, 4)
ax3.hist2d(-Z, Y, bins=(50,50))
ax3.set_xlabel('z (cm)')
ax3.set_ylabel('y (cm)')

plt.show()

^ Minimum working example ^

I have tried the surface_plot method, but none of my data is in the form of a 2D array for the z argument.

Here is an image I am currently producing, it is hard to produce a minimal example from the database I am working with, but basically I want to format all of these subplots onto one plot.

You can calculate the information of ax.hist2d via np.histogram2d and then use these values to color a mesh which you can position in 3D.

Here is some code to demonstrate the idea. This code also inverts the x-axis and uses the order X, -Z, Y for the scatter plot, as in the original question. As the plot looks quite busy, the scatter plot might be left out.

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(20230223)
M, N = 2000, 6
X = np.random.normal(0, 2, size=(M, N)).cumsum(axis=0).ravel() + 5
Y = np.random.normal(0, 2, size=(M, N)).cumsum(axis=0).ravel() + 5
Z = np.random.normal(0, 2, size=(M, N)).cumsum(axis=0).ravel() + 5

fig = plt.figure(figsize=(12, 12))

ax = fig.add_subplot(1, 1, 1, projection='3d')
ax.scatter(X, -Z, Y, alpha=0.9, s=2)
ax.invert_xaxis()
ax.set_xlabel('x (cm)')
ax.set_ylabel('z (cm)')
ax.set_zlabel('y (cm)')

cmap = plt.cm.YlOrRd

counts, xbins, ybins = np.histogram2d(-Z, Y, bins=(50, 50))
yg, xg = np.meshgrid(ybins, xbins)
norm = plt.Normalize(vmin=1, vmax=counts.max())
counts[counts == 0] = np.nan
ax.plot_surface(np.full_like(xg, X.max()), xg, yg, facecolors=cmap(norm(counts)), rstride=1, cstride=1, shade=False)

counts, xbins, ybins = np.histogram2d(X, Y, bins=(50, 50))
yg, xg = np.meshgrid(ybins, xbins)
norm = plt.Normalize(vmin=1, vmax=counts.max())
counts[counts == 0] = np.nan
ax.plot_surface(xg, np.full_like(xg, -Z.min()), yg, facecolors=cmap(norm(counts)), rstride=1, cstride=1, shade=False)

counts, xbins, ybins = np.histogram2d(X, -Z, bins=(50, 50))
yg, xg = np.meshgrid(ybins, xbins)  # bins as a 2D grid
norm = plt.Normalize(vmin=1, vmax=counts.max())
counts[counts == 0] = np.nan  # make 0 transparent
ax.plot_surface(xg, yg, np.full_like(xg, Y.min()), facecolors=cmap(norm(counts)), rstride=1, cstride=1, shade=False)

plt.show()