I have a discrete 3d mass density distribution function, $\rho(x, y, z)$, representing the distribution of mass in a sphere with known radius $R_{{0}}$ centered at ($0$, $0$). In other words, given any 3d position inside this sphere ($x$, $y$, $z$), taken to be the center of some cube grid with a fixed otherwise small size d, this 3-variable function would give me the mass density inside that grid.
The mass happens to be concentrated in few clumps across the sphere the size of which are much larger than the size of the cube grid.
This 3d system is projected in xy plane as seen from +z axis. Hence, the densest regions of this sphere is also translated into a 2d region along the axis of projection. I happen to have the array of 2d positions of the square grids corresponding to this 2d dense region.
Now, I would like to utilized these two pieces of information, namely 3d density distribution function and 2d array of the densest region in one projection, to calculate the 1d distribution of mass density along the beam passing through any one of the grids of the 2d array perpendicular to the xy plane along the third dimension. Once I have that 1d mass density distribution function at any point in xy plane, I can find the position of the median of the distribution. Since my function is discrete, I am not sure how to approach it numerically. But if I get some help in the analysis of the problem, I should be able to write a code in python to do the rest.
To go from a joint probability distribution function to a marginal distribution function (what you seem to be describing) then you integrate along the dimensions you want to eliminate. In this case, since you only have points on a grid, you can approximate that integration as a Riemann sum. You'll probably want to normalize afterwards to make sure that your end product integrates to 1 (again, you'll probably have to approximate this integration as a Riemann sum).