matlab 抽样数量,基于matlab的多元定制累积分布函数抽样

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