Generate Y from the conditional f(y|x)

I want to sample from a 3 dimensional probability density function. The method that I am using is to discretize the density funcion by evaluating it at regular intervals. After normalizing, I have a three dimensional matrix of the "probability" of finding a particle at a given point along x, y, z coordinates. Then I calculated the marginal of x by summing the probabilities along y and z. I use the method of inverse transform sampling to generate n x coordinates from the marginal of x.

Now I need to generate n y coordinates GIVEN the x coordinates (of course each of the n points will have a different x coordinate). I do not know of an efficient way to do this; I tried summing the probabilities of y coordinates along the z coordinates to get the probability of a y coordinate as a function of a given x coordinate (I don't think we can call this the marginal of Y, but it is similar). Now I tried to use use 'cumsum' to get the cumulative probability of y coordinates along each x coordinate. The problem is that now each sample must be taken from a different CDF (because the probability of distribution of y coordinates depends upon the already given x coordinates).

[Peter]

Here is my code:

b = .6;
  width = 1*10^-6;
  height = 1*10^-6;
  numpart = 10;
  %-------------------
  % Generate X from the marginal f(x), Y from the conditional f(y|x), and Z from the conditional f(z|x & y)"
  PDF = @(x,y,z) 1 - exp( -b*exp(-4*x.^2/width^2 -4*y.^2/width^2 -4*z.^2/height^2 ) );  % the distribution of points

% discretize the PDF
x = -1*width:width/2:1*width;
y = -1*width:width/2:1*width;
z = -1*height:height/2:1*height;

[x,y,z] = meshgrid(x,y,z);
PDFdiscrete = PDF(x,y,z);
norm = sum(PDFdiscrete(:));    % normalize so that the sums of the marginals are unity
PDFdiscrete = PDFdiscrete/norm;  % the discrete PDF

% calculate the marginal of X
MargLength = size(PDFdiscrete,2);
PDFperm = ipermute(PDFdiscrete, [1 3 2]);
m = reshape(PDFperm, [], MargLength);
MarginalX = sum(m, 1);  % sum along y and z for each value of x
% sum(MarginalX) % check that this equals 1
[~ , draw_samples] = histc(rand(numpart , 1), [0 cumsum(MarginalX)]);
position(:,1) = draw_samples .*sizeX -width;

% calculate the marginal of Y as a function of X
MarginalY = sum(PDFdiscrete, 3);  % sum along y and z for each value of x
cumMargY = cumsum(MarginalY);
[~ , draw_samplesY] = histc(rand(numpart , 1), [0 cumMargY(:,draw_samples)']);  % my problem is at this point: the prob distribution is different for each of the sampled points
position(:,2) = draw_samplesY .*sizeY -width;

Answer 1

[Teja Muppirala]

If your goal is to generate points with that 3-dimensional PDF, then I think it could be done a bit simpler without having to do all sorts of cumbersome manipulations involving marginal distributions. Generate uniform random numbers, and then remove the ones that don't fit "under" the probability distribution.

In other words, a generalization of Richard Brown's answer here: http://www.mathworks.com/matlabcentral/answers/35861-how-to-sample-from-custom-2d-distribution

b = .6;
width = 1*10^-6;
height = 1*10^-6;

P = @(x,y,z) 1 - exp( -b*exp(-4*x.^2/width^2 -4*y.^2/width^2 -4*z.^2/height^2 ) );  % the distribution of points
Pmax = P(0,0,0); % Maximum value in the PDF

% Guess some x,y,z values
Ncandidates = 1e6;
X = ( -3+6.*rand(Ncandidates,1) )*width;
Y = ( -3+6.*rand(Ncandidates,1) )*width;
Z = ( -3+6.*rand(Ncandidates,1) )*height;
Pguess = Pmax*rand(Ncandidates,1); %Guess a probability value

good = find( Pguess < P(X,Y,Z) ); %Find the ones "under" the PDF

plot3(X(good),Y(good),Z(good),'.');
axis(1e-6*[-2 2 -2 2 -2 2])
axis vis3d
title([num2str(numel(good)) ' Points with the given distribution']);
box on

In this particular case, you could also take advantage of the fact that it is radially symmetric (R^2 = X^2 + Y^2 + Z^2), so you would only need to generate radii, and then set the direction randomly over a sphere:

http://www.mathworks.com/matlabcentral/answers/44481#answer_54528

[Peter]

This works, but I don't fully understand why. Why are the random candidates multiplied by 6 and subtracted by 3?