Clear Filters
Clear Filters

Generate Y from the conditional f(y|x)

4 views (last 30 days)
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).
  1 Comment
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;

Sign in to comment.

Accepted Answer

Teja Muppirala
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.
In other words, a generalization of Richard Brown's answer here:
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
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:
  1 Comment
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?

Sign in to comment.

More Answers (0)


Find more on Random Number Generation in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!