# Generate Y from the conditional f(y|x)

7 views (last 30 days)
Peter on 17 Jul 2013
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 on 17 Jul 2013
Edited: Peter on 17 Jul 2013
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;

Teja Muppirala on 18 Jul 2013
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.
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:
Peter on 18 Jul 2013
Edited: Peter on 18 Jul 2013
This works, but I don't fully understand why. Why are the random candidates multiplied by 6 and subtracted by 3?