Generate Y from the conditional f(y|x)

3 次查看（过去 30 天）
Peter
Peter 2013-7-17
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 个评论
Peter
Peter 2013-7-17
编辑：Peter 2013-7-17
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
Teja Muppirala 2013-7-18
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
  1 个评论
Peter
Peter 2013-7-18
编辑：Peter 2013-7-18
This works, but I don't fully understand why. Why are the random candidates multiplied by 6 and subtracted by 3?

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Probability Distributions 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by