从 Matlab 中的多元自定义累积分布函数中抽样
Sampling from multivariate customised cumulative distribution function in Matlab
考虑一个 5 变量 cumulative distribution function (cdf)，我称之为 F。
我想在 Matlab 中从这个 cdf 中随机采样 5x1 向量。 F不是已经在Matlab中实现过的cdf（比如normal，t-student等）。具体定义为
我在这个论坛和其他论坛上阅读了一些关于如何从 Matlab 中的自定义概率分布函数中采样的文章 question/answers。然而，
1) 其中大部分是针对单变量 cdf 的，例如. The idea is to apply inverse transform sampling。我的问题有点复杂，因为我需要 "invert" 一个 5 变量函数。
2) 另一种选择是使用 slicesample as suggested here 但我不知道如何在我的案例中编写概率密度函数的解析表达式。
3) Here 的另一个想法，但特定于双变量情况。
能否请您帮助我了解如何继续？
Your link under #3 给出了解决方案的提示。它解释了当您拥有 PDF 时的双变量情况。在这里，我们将把它扩展到任意数量的维度，适用于你有 CDF 的情况。
所以过程是：
- 计算 r1.
的边缘 CDF
- 使用这个边际 CDF 随机抽样（
解释了如何做到这一点）。
- 计算 r2 给定 r1 的边缘 CDF .
- 使用这个边际 CDF 随机抽样
- 计算 r3 给定 r1 的边缘 CDF 和 r2。
- 等等等。你知道这是怎么回事。
请注意，如果您有 PDF，计算边际分布涉及对剩余变量进行积分。因此 r1 的边际分布需要对 r2[=151 进行积分=]..r5，r2[=的边际分布151=] 给定 r1 需要在 r3[=151= 上积分..r5等
当你有 CDF 时，计算边际分布很简单，因为它已经集成了 PDF：r1[=151= 的边际分布] 是 F(x,∞,∞,∞,∞)。但是，在给定一个或多个变量的情况下获得边际分布需要微分：给定 r[=150]r2 的边际分布=]1 需要沿 r1 微分，r[ 的边际分布=150=]3 给定 r1 和 r 2 需要区分 r1 和 r2,等等
也许可以通过分析获得这些导数（这将是更有效的解决方案）。这里我们使用有限差分导数近似代替（这使得插入任何 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 选择了一些随机值，并使用它们定义了一个包含 5 个变量 r1..[=20= 的函数 F ].我认为 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
这包含相当多的复杂性。它在给定固定值 r(1:ii-1).
的点 x 沿给定维度 ii 对 CDF 进行采样
复杂性来自计算偏导数。如果我们要在没有选择任何固定值的情况下计算任何一维的边际分布，我们会简单地做例如
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 是一个包含函数 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
考虑一个 5 变量 cumulative distribution function (cdf)，我称之为 F。
我想在 Matlab 中从这个 cdf 中随机采样 5x1 向量。 F不是已经在Matlab中实现过的cdf（比如normal，t-student等）。具体定义为
我在这个论坛和其他论坛上阅读了一些关于如何从 Matlab 中的自定义概率分布函数中采样的文章 question/answers。然而，
1) 其中大部分是针对单变量 cdf 的，例如
2) 另一种选择是使用 slicesample as suggested here 但我不知道如何在我的案例中编写概率密度函数的解析表达式。
3) Here 的另一个想法，但特定于双变量情况。
能否请您帮助我了解如何继续？
Your link under #3 给出了解决方案的提示。它解释了当您拥有 PDF 时的双变量情况。在这里，我们将把它扩展到任意数量的维度，适用于你有 CDF 的情况。
所以过程是：
- 计算 r1. 的边缘 CDF
- 使用这个边际 CDF 随机抽样（
- 计算 r2 给定 r1 的边缘 CDF .
- 使用这个边际 CDF 随机抽样
- 计算 r3 给定 r1 的边缘 CDF 和 r2。
- 等等等。你知道这是怎么回事。
请注意，如果您有 PDF，计算边际分布涉及对剩余变量进行积分。因此 r1 的边际分布需要对 r2[=151 进行积分=]..r5，r2[=的边际分布151=] 给定 r1 需要在 r3[=151= 上积分..r5等
当你有 CDF 时，计算边际分布很简单，因为它已经集成了 PDF：r1[=151= 的边际分布] 是 F(x,∞,∞,∞,∞)。但是，在给定一个或多个变量的情况下获得边际分布需要微分：给定 r[=150]r2 的边际分布=]1 需要沿 r1 微分，r[ 的边际分布=150=]3 给定 r1 和 r 2 需要区分 r1 和 r2,等等
也许可以通过分析获得这些导数（这将是更有效的解决方案）。这里我们使用有限差分导数近似代替（这使得插入任何 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 选择了一些随机值，并使用它们定义了一个包含 5 个变量 r1..[=20= 的函数 F ].我认为 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
这包含相当多的复杂性。它在给定固定值 r(1:ii-1).
x 沿给定维度 ii 对 CDF 进行采样
复杂性来自计算偏导数。如果我们要在没有选择任何固定值的情况下计算任何一维的边际分布，我们会简单地做例如
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 是一个包含函数 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