Generate Y from the conditional f(y|x)
16 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).
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: 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
axis(1e-6*[-2 2 -2 2 -2 2])
title([num2str(numel(good)) ' Points with the given distribution']);
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: