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

Question

考虑一个我称之为 F 的 5 变量 cumulative distribution function (cdf)。

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

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

1) 其中大多数是用于单变量 cdf，例如here。这个想法是申请inverse transform sampling。我的问题有点复杂，因为我需要“反转”一个 5 变量函数。

2) 另一种选择是使用slicesample 建议here，但我不知道如何编写概率密度函数的解析表达式。

3) Here 的另一个想法，但特定于双变量情况。

您能否帮助我了解如何进行？

Answer 1

Your link under #3 给出了解决方案的提示。它解释了当您拥有 PDF 时的双变量情况。在这里，我们会将其扩展到任意数量的维度，以应对您拥有 CDF 的情况。

所以流程是：

1. 计算 r₁ 的边际 CDF。
2. 使用此边缘 CDF 进行随机抽样（another link you posted 解释了如何执行此操作）。
3. 在给定 r₁ 的情况下计算 r₂ 的边际 CDF。
4. 使用此边缘 CDF 随机抽样
5. 在给定 r₁ 和 r 的情况下计算 r₃ 的边际 CDF ₂。
6. 等。等等。你知道这是怎么回事。

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

当你有一个 CDF 时，计算边际分布是微不足道的，因为它已经集成了 PDF：r₁ 的边际分布是 F em>(x,∞,∞,∞,∞)。但是，在给定一个或多个变量的情况下获得边际分布需要区分： r₂ 的边际分布给定 r₁ 需要沿 r₁ 进行微分，给定 r 的 r₃ 的边际分布₁ 和 r₂ 需要区分 r₁ 和 r₂ 等

也许可以通过分析获得这些导数（这将是更有效的解决方案）。这里我们改用有限差分导数逼近（这样更容易插入任何 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 和sigma_b 选择了一些随机值，并使用它们定义了一个函数F 的5 个变量r1..r5。我认为 PDF 的域在所有维度上都是相同的，我发现一个区域比实际需要的稍大（lims）。接下来，我通过调用sample_random从分布F中获取1000个随机样本：

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 的最大值之间选择一个随机值（请注意，当我们为每个维度选择 r 的值时，最大值会减小，当 ii==1 最大值为 1 时），我们使用这个随机值来插值到逆采样的边际 CDF 中（逆简单地意味着交换 x 和 y）。我需要从marginal 中删除重复值，因为它变成了interp1 中的x，并且此函数要求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 固定值执行此操作。这需要对这些固定值中的每一个以及delta 移位的每个组合采样F 两次，总共n^2 次（对于n 固定值）。 dec2bin 位用于获取所有这些组合。 sign 确定是从运行总数中添加还是减去给定样本。 args 是一个元胞数组，函数 F 有 5 个参数，元素 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