时间戳

首页 > 解决方案 > 如何修复斜梯形分布采样输出样本大小

algorithm - 如何修复斜梯形分布采样输出样本大小

问题描述

我正在尝试使用逆变换采样生成倾斜的梯形分布。输入是斜坡开始和结束的值（a、b、c、d）和样本大小。

a=-3;b=-1;c=1;d=8; 
SampleSize=10e4;
h=2/(d+c-a-b);

然后我计算斜坡和平面组件的长度之比，以获得每个的样本大小：

firstramp=round(((b-a)/(d-a)),3);
flat=round((c-b)/(d-a),3);
secondramp=round((d-c)/(d-a),3);
n1=firstramp*SampleSize; %sample size for first ramp
n3=secondramp*SampleSize; %sample size for second ramp
n2=flat*SampleSize;

最后我从以下代码中得到直方图：

quartile1=h/2*(b-a);
quartile2=1-h/2*(d-c);

y1=linspace(0,quartile1,n1);
y2=linspace(quartile1,quartile2,n2);
y3=linspace(quartile2,1,n3);

%inverse cumulative distribution functions
invcdf1=a+sqrt(2*(b-a)/h)*sqrt(y1);
invcdf2=(a+b)/2+y2/h;
invcdf3=d-sqrt(2*(d-c)/h)*sqrt(1-y3);

distr=[invcdf1 invcdf2 invcdf3];

histogram(distr,100)

然而，斜坡和平坦分量的采样不相等，如下所示：

混乱的输出

我通过反复试验解决了这个问题，将坡道的样本量减少了一半：

n1=0.5*firstramp*SampleSize; %sample size for first ramp
n3=0.5*secondramp*SampleSize; %sample size for second ramp
n2=flat*SampleSize;

这使得分布看起来像这样：

固定分布

然而，这使得输出样本小于输入中给出的样本。

我还尝试了改变坡道和平坦的样本大小的不同组合。这也有效：

n1=0.75*firstramp*SampleSize; %sample size for first ramp
n3=0.75*secondramp*SampleSize; %sample size for second ramp
n2=1.5*flat*SampleSize;

它增加了输出样本，但仍然没有接近。

任何帮助将不胜感激。

完整代码：

a=-3;b=-1;c=1;d=8; 
SampleSize=10e4;%*1.33333333333333;
h=2/(d+c-a-b);
firstramp=round(((b-a)/(d-a)),3);
flat=round((c-b)/(d-a),3);
secondramp=round((d-c)/(d-a),3);

n1=firstramp*SampleSize; %sample size for first ramp
n3=secondramp*SampleSize; %sample size for second ramp
n2=flat*SampleSize;

quartile1=h/2*(b-a);
quartile2=1-h/2*(d-c);

y1=linspace(0,quartile1,.75*n1);
y2=linspace(quartile1,quartile2,1.5*n2);
y3=linspace(quartile2,1,.75*n3);

%inverse cumulative distribution functions
invcdf1=a+sqrt(2*(b-a)/h)*sqrt(y1);
invcdf2=(a+b)/2+y2/h;
invcdf3=d-sqrt(2*(d-c)/h)*sqrt(1-y3);

distr=[invcdf1 invcdf2 invcdf3];

histogram(distr,100)
%end

标签： algorithmmatlabmathdistributionskew

解决方案

我不知道 Matlab，所以我希望其他人能参与进来，但由于这里没有人这样做。

如果我正确阅读了您的代码，那么您所做的并不是倒置。反转是 1-1，即一个统一的输入产生一个结果。您似乎正在使用一种称为“组合方法”的技术。在组合中，整体分布由组件组成，每个组件都可以直接生成。您可以根据它们相对于整体的比例/概率来选择从哪个组件生成。对于密度函数，概率是密度曲线下的面积，因此您的第一个错误是相对于每个分量的宽度而不是使用它们的面积来采样分量。firstramp对于您指定的、flat和，正确的抽样比例是 2/13、4/13 和 7/13secondramp组件，分别。第二个错误（相对较小）是将准确的样本大小分配给每个组件。概率为 2/13 并不意味着2*SampleSize/13您的样本完全来自firstramp，这意味着这是该组件的预期样本大小。随机变量的期望值不一定（甚至可能是）您实际得到的结果。

在伪代码中，组合方法是

generate U ~ Uniform(0,1)
if U <= 2/13:
   generate and return a value from firstramp
else if U <= 6/13:
   generate and return a value from flat
else:
   generate and return a value from secondramp

请注意，由于每个generate选项都将使用一个或多个制服，并且在选项之间进行选择需要制服U，因此这不是倒置。

如果你想要一个实际的反演，你需要量化你的密度，整合它以获得累积分布函数，然后通过设置F(X) = U和求解来应用反演技术X。由于您的分布由不同的组件组成，因此密度和累积密度都是分段函数。

在根据两个三角形的面积和平面部分的面积必须加起来为 1 的要求推导出高度后，我为您的密度得出以下结论：

       | (x + 3) / 13       -3 <= x <= -1
       |
f(x) = | 2 / 13             -1 <= x <= 1
       |
       | 2 * (8 - x) / 91    1 <= x <= 8

将这个和收集项整合起来会产生 CDF：

       | (x + 3)**2 / 26                    -3 <= x <= -1
       |
F(x) = | (2 + x) * 2 / 13                   -1 <= x <= 1
       |
       | 6 / 13 + [49 - (x - 8)**2] / 91     1 <= x <= 8

最后，确定段之间断点处的 F(x) 值并应用反演产生以下伪代码算法：

generate U ~ Uniform(0,1)
if U <= 2 / 13:
    return 2 * sqrt( (13 * U) / 2 ) - 3
else if U <= 6 / 13:
    return (13 * U) / 2 - 2:
else:
    return 8 - sqrt( 91 * (1 - U) )

请注意，这是一个真正的反转。结果是通过生成单个 U 来确定的，并根据它所在的范围以不同的方式对其进行转换。

推荐阅读