Your link under #3
给出解决方案的提示。它解释了当你有一个pdf文件时的二元情况。在这里,我们将把这个扩展到任意数量的维度,当你有一个CDF的时候。
所以这个过程是:
计算边际cdf
R
一
.
使用这个边缘cdf随机抽取一个样本(
another link you posted
说明如何执行此操作)。
计算边际cdf
R
二
鉴于
R
一
.
使用这个边缘cdf随机抽取一个样本
计算边际cdf
R
三
鉴于
R
一
和
R
二
.
等等等等。你知道这是怎么回事了。
请注意,如果你有一个PDF,计算边际分布涉及整合剩余的变量。因此,边际分布
R
一
需要通过
R
二
…
R
五
,边际分布
R
二
鉴于
R
一
需要通过
R
三
…
R
五
等。
当您有一个cdf时,计算边际分布是很简单的,因为它已经集成了pdf:
R
一
是
f
(
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
)接下来,我从分布中获得1000个随机样本
f
通过呼叫
sample_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的最大值之间的随机值(注意,最大值随着我们选择
R
对于每个维度,当
ii==1
最大值是1),我们使用这个随机值插值到逆采样边缘cdf(逆简单地说就是交换x和y)。我需要从
marginal
因为它变成了
X
在里面
interp1
,此函数需要
X
值是唯一的。
最后,函数
小精灵
:
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
这包含了相当多的复杂性。它沿着给定的维度对cdf进行采样
二
,在点
X
,给定固定值
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);
(这是对第一维偏导数的近似)。
代码在
小精灵
这是为了
ii-1
固定值。这需要取样
f
对于每个固定值和每个组合
三角洲
轮班,总共
n^2
时代(对)
n
固定值)。这个
dec2bin
比特是得到所有这些组合。
sign
确定是否从正在运行的总数中添加或减去给定的样本。
args
是一个具有5个函数参数的单元格数组
f
元素
1:ii-1
是固定值,元素
二
设置为
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