从 Matlab 中的多元自定义累积分布函数中采样

Question

考虑一个 5 变量累积分布函数(cdf)，我将其称为 F。

我想在 Matlab 中从这个 cdf 中随机采样 5x1 向量。F不是Matlab中已经实现的cdf（如normal、t-student等）。具体来说，它被定义为

我已经在这个论坛和其他论坛上阅读了几个关于如何从 Matlab 中自定义概率分布函数进行采样的问题/答案。然而，

1）大多数都是针对单变量cdf，例如这里。这个想法是应用逆变换采样。我的问题有点复杂，因为我需要“反转”5 变量函数。

2）另一种选择可能是按照此处的建议使用slicesample，但我不知道在我的情况下如何编写概率密度函数的解析表达式。

3）这是另一个想法，但特定于二元情况。

您能帮助我了解如何继续吗？

Answer 1

#3 下的链接给出了解决方案的提示。它解释了当您拥有 PDF 时的双变量情况。在这里，我们将把它扩展到任意数量的维度，以防您有 CDF 的情况。

所以流程是：

1. 计算r ₁的边际 CDF 。
2. 使用此边际 CDF 进行随机采样（您发布的另一个链接解释了如何执行此操作）。
3. 计算给定r _1的r ₂的边际 CDF 。
4. 使用边际 CDF 进行随机样本
5. 给定r ₁和r _{2 ，计算}r ₃的边际 CDF 。
6. 等等等等。你知道这是怎么回事了。

请注意，如果您有 PDF，计算边际分布涉及对剩余变量进行积分。因此， r ₁的边际分布需要对r ₂ .. r ₅进行积分，给定r _1的r _{2 的}边际分布需要对r ₃ .. r ₅进行积分，等等。

当您拥有 CDF 时，计算边际分布是微不足道的，因为它已经对 PDF 进行了积分：r ₁的边际分布为F ( x , ∞ , ∞ , ∞ , ∞ ) 。然而，获得给定一个或多个变量的边际分布需要微分：给定r _1的r ₂边际分布需要沿r ₁微分，给定r ₁和r _2的r ₃边际分布需要沿r ₁和r ₂微分， ETC。

或许可以通过分析获得这些导数（这将是更有效的解决方案）。这里我们使用有限差分导数近似（这使得插入任何 CDF 变得更容易）。

让我们看一些 MATLAB 代码：

sigma_a = 0.5;
sigma_b = 0.3;
F = @(r1,r2,r3,r4,r5)exp(-exp(-r1) - (exp(-r2/sigma_a)+exp(-r3/sigma_a)).^sigma_a ...
                                   - (exp(-r4/sigma_b)+exp(-r5/sigma_b)).^sigma_b);

lims = [-5,10]; % This is the area along all dimensions containing 99.99% of the PDF

N = 1000;
values = zeros(N,5);
for n=1:N
   values(n,:) = sample_random(F,5,lims);
end

在这里，我为sigma_a和选取了一些随机值，并使用这些值定义了5 个变量的sigma_b函数...... 我认为 PDF 的域在所有维度上都是相同的，我发现一个区域比实际需要的稍大 ( )。接下来，我通过调用以下命令从分布中获取 1000 个随机样本：Fr1r5limsFsample_random

function r = sample_random(F,N,lims)
delta = diff(lims)/10000;
x = linspace(lims(1),lims(2),300);
r = inf(1,N);
for ii = 1:N
   marginal = get_marginal(F,r,ii,x,delta);
   p = rand * marginal(end);
   [~,I] = unique(marginal); % interp1 cannot handle duplicated points, let's remove them
   r(ii) = interp1(marginal(I),x(I),p);
end

delta是我们用于导数的有限差分近似的距离。x表示沿任意一维的样本点F。

我们首先定义r为向量[inf,inf,inf,inf,inf]，我们将使用它作为样本位置，并且在函数末尾它将包含从我们的分布中提取的随机值。

接下来，我们循环 5 个维度，在每次迭代中，我们对维度 的边际分布进行采样ii，并给出先前维度的值（已被选取）。该函数get_marginal如下。我们在 0 和该边际 CDF 的最大值之间选择一个随机值（请注意，当最大值为 1r时，当我们为每个维度选择 的值时，最大值会减小ii==1），并且我们使用该随机值插值到逆采样边际中CDF（逆简单来说就是交换 x 和 y）。我需要删除重复的值 frommarginal因为它变成了xin interp1，并且此函数要求x值是唯一的。

最后，函数get_marginal：

function marginal = get_marginal(F,r,ii,x,delta)
N = length(r);
marginal = zeros(size(x));
for jj=0:2^(ii-1)-1
   rr = flip(dec2bin(jj,N)-'0');
   sign = mod(sum(rr,2),2);
   if sign == 0
      sign = 1;
   else
      sign = -1;
   end
   args = num2cell(r - delta * rr);
   args{ii} = x;
   marginal = marginal + sign * F(args{:});
end

这包含相当多的复杂性。ii它在给定固定值的点上沿x给定维度对 CDF 进行采样r(1:ii-1)。

复杂性来自于计算偏导数。如果我们要在尚未选择任何固定值的情况下计算任何一维的边际分布，我们只需执行以下操作：

marginal = F(inf,x,inf,inf,inf);

选择一个值后，我们会做

marginal = F(r1,x,inf,inf,inf) - F(r1-delta,x,inf,inf,inf);

（这是沿第一维的偏导数的近似值）。

中的代码get_marginal对ii-1固定值执行此操作。F这需要对每个固定值以及每个delta班次组合采样两次n^2（对于n固定值）。关键dec2bin是获得所有这些组合。sign确定是否从运行总计中添加或减去给定样本。args是一个元胞数组，其中包含函数 的 5 个参数F，元素1:ii-1是固定值，元素ii设置为x，元素ii+1:N为inf。

最后，我绘制数据集的边际分布values（其中包含从 CDF 中随机抽取的 1000 个元素），并与 CDF 的边际分布叠加，以验证所有内容是否正确：

lims = [-2,5];
x = linspace(lims(1),lims(2),300);
figure
for ii=1:5
   subplot(5,1,ii)
   histogram(values(:,ii),'normalization','cdf','BinLimits',lims)
   hold on
   args = num2cell(inf(1,5));
   args{ii} = x;
   plot(x,F(args{:}))
   text(5.2,0.5,['r_',num2str(ii)])
end