在Matlab中从任意离散概率密度函数生成随机样本

Question

我有一个在Matlab中以矩阵形式离散化的任意概率密度函数，这意味着对于每个x,y对，概率都存储在矩阵中： A(x,y) = 概率

这是一个100x100的矩阵，我想能够从该矩阵生成二维随机样本(x,y)，并且如果可能的话，能够计算PDF的均值和其他矩。我想这样做是因为在重新采样后，我想将样本拟合到近似的高斯混合模型。

我已经到处查找，但没有找到像这样具体的内容。希望您能帮助我。

谢谢。

- Miguel

---

我不能给你代码。但是如果在文档中找不到某些东西，你可以自己实现它。你只需要能够从一个离散分布中进行采样。这篇Wiki-Article展示了一些方法，其中有一些非常容易实现！如果速度不那么重要：选择线性搜索。如果速度很重要：选择Alias-Method。 - sascha

我认为这个问题不应该在这里提问。从任意概率密度函数中计算均值和其他矩总是很困难的，但如果你能获得条件概率：x|y和y|x，那么你就可以使用“吉布斯采样”来得到你想要的结果。你可以在这里找到一个例子（http://timsalimans.com/the-power-of-jit-compilation/）。 - Gnimuc

Answer 1

如果您真的有一个由A定义的离散概率密度函数（而不是仅由A描述的连续概率密度函数），那么您可以通过将二维问题转换为一维问题来“作弊”。

%define the possible values for the (x,y) pair
row_vals = [1:size(A,1)]'*ones(1,size(A,2));  %all x values
col_vals = ones(size(A,1),1)*[1:size(A,2)];  %all y values

%convert your 2D problem into a 1D problem
A = A(:);
row_vals = row_vals(:);
col_vals = col_vals(:);

%calculate your fake 1D CDF, assumes sum(A(:))==1
CDF = cumsum(A); %remember, first term out of of cumsum is not zero

%because of the operation we're doing below (interp1 followed by ceil)
%we need the CDF to start at zero
CDF = [0; CDF(:)];

%generate random values
N_vals = 1000;  %give me 1000 values
rand_vals = rand(N_vals,1);  %spans zero to one

%look into CDF to see which index the rand val corresponds to
out_val = interp1(CDF,[0:1/(length(CDF)-1):1],rand_vals); %spans zero to one
ind = ceil(out_val*length(A));

%using the inds, you can lookup each pair of values
xy_values = [row_vals(ind) col_vals(ind)];

我希望这能帮到你！ 芯片

- chipaudette

Answer 2

我不认为Matlab具备生成任意分布多元随机变量的内置功能。事实上，对于一元随机数同样如此。但是，虽然可以基于累积分布函数轻松生成后者，但多元分布不存在CDF，因此生成这些数字要困难得多（主要问题在于2个或更多变量之间存在相关性）。因此，你提出的这个问题超出了本网站的范围。

既然半个答案胜过没有答案，以下是使用Matlab进行数值计算的平均值和更高次矩的方法：

%generate some dummy input
xv=linspace(-50,50,101);
yv=linspace(-30,30,100);
[x y]=meshgrid(xv,yv);

%define a discretized two-hump Gaussian distribution
A=floor(15*exp(-((x-10).^2+y.^2)/100)+15*exp(-((x+25).^2+y.^2)/100));
A=A/sum(A(:)); %normalized to sum to 1

%plot it if you like
%figure;
%surf(x,y,A)

%actual half-answer starts here    

%get normalized pdf
weight=trapz(xv,trapz(yv,A));
A=A/weight; %A normalized to 1 according to trapz^2

%mean
mean_x=trapz(xv,trapz(yv,A.*x));
mean_y=trapz(xv,trapz(yv,A.*y));

因此，关键是您可以使用两个连续的trapz调用在矩形网格上执行双重积分。这使您可以计算任何具有与网格相同形状但向量组件必须独立计算的量的积分。如果您只希望计算可以用x和y参数化的内容（这些内容自然与您的网格大小相同），则可以不进行任何额外思考而完成。您还可以为集成定义一个函数：

function res=trapz2(xv,yv,A,arg)

if ~isscalar(arg) && any(size(arg)~=size(A))
    error('Size of A and var must be the same!')
end

res=trapz(xv,trapz(yv,A.*arg));

end

这样您就可以计算诸如以下内容：

weight=trapz2(xv,yv,A,1);
mean_x=trapz2(xv,yv,A,x);

注意：我在示例中使用101x100网格的原因是双重调用trapz应按正确顺序执行。如果在调用中交换xv和yv，则由于与A的定义不一致，您将得到错误的答案，但如果A是正方形，则这不会明显。我建议在开发阶段避免对称量。

- Andras Deak -- Слава Україні

---

选择101x100网格而不是100x100网格是个好主意。在整个代码中正确地获取尺寸和形状非常棘手（但非常重要）。使数组不是正方形是一个很好的方法来做到这一点！ - chipaudette