1

我想计算两个三角随机变量的总和，

P(x1+x2 < y)

图片

有没有更快的方法在 Matlab 中实现两个三角形随机变量的总和？

编辑：似乎有一种更简单的方法，如这个minitab演示中所示。所以这不是不可能的。遗憾的是，它没有解释如何计算 PDF。仍在研究如何在 matlab 中做到这一点。

EDIT2：根据建议，我正在使用convMatlab 中的函数来开发两个随机变量之和的 PDF：

clear all;
clc;


pd1 = makedist('Triangular','a',85,'b',90,'c',100);
pd2 = makedist('Triangular','a',90,'b',100,'c',110);

x = linspace(85,290,200);
x1 = linspace(85,100,200);
x2 = linspace(90,110,200);
pdf1 = pdf(pd1,x1);
pdf2 = pdf(pd2,x2);

z = median(diff(x))*conv(pdf1,pdf2,'same');

p1  = trapz(x1,pdf1) %probability P(x1<y)
p2  = trapz(x2,pdf2) %probability P(x2<y)
p12 = trapz(x,z)     %probability P(x1+x2 <y)

hold on;
plot(x1,pdf1) %plot pdf of dist. x1
plot(x2,pdf2) %plot pdf of dist. x2
plot(x,z)     %plot pdf of x1+x2
hold off;

但是这段代码有两个问题：

  1. X1+X2 的 PDF 积分远高于 1。
  2. X1+X2 的 PDF 根据 x 的范围变化很​​大。直观地说，如果 X1+X2 大于 210（两个独立三角形分布的上限“c”之和，100 + 110），那么 P(X1+X2 <210) 不应该等于 1 吗？此外，由于下限“a”是 85 和 90，P(X1+X2 <85) = 0？
4

2 回答 2

2

自变量之和的 pdf 是变量 pdf 的卷积。所以你需要用三角pdf计算两个变量的卷积。三角形是分段线性的，因此卷积将是分段二次的。

有几种方法可以解决它。如果数值结果是可接受的：离散化 pdf 并计算离散化 pdf 的卷积。我相信convMatlab中有一个功能。如果没有，您可以进行快速傅里叶变换（通过fft），逐点计算乘积，然后进行逆变换（ifft如果我没记错的话），因为 fft(convolution(f, g)) = fft(f) fft(g ）。conv如果您使用或，则需要小心规范化fft。

如果你必须有一个精确的结果，卷积只是一个积分，如果你小心积分的限制，你可以手动计算出来。我不知道您是否可以使用 Matlab 符号工具箱，如果可以，我不知道它是否可以处理分段定义的函数的积分。

于 2015-12-21T20:30:02.270 回答
1

以下是未来用户的正确实施。非常感谢 Robert Dodier 的指导。

clear all;
clc;

min1 = 85;
max1 = 100;
min2 = 90;
max2 = 110;
y    = 210;

pd1 = makedist('Triangular','a',min1,'b',90,'c',max1);
pd2 = makedist('Triangular','a',min2,'b',100,'c',max2);

dx = 0.01; % to ensure constant spacing
x1 = min1:dx:max1; % Could include some of the region where
x2 = min2:dx:max2; % the pdf is 0, but we don't have to.

x12 = linspace(...
    x1(1)   + x2(1) , ...
    x1(end) + x2(end) , ...
    length(x1)+length(x2)-1);
[c,index] = min(abs(x12-y));
x_short = linspace(min1+min2,x12(index),index);

pdf1 = pdf(pd1,x1);
pdf2 = pdf(pd2,x2);
pdf12 = conv(pdf1,pdf2)*dx;

zz = pdf12(1:index);
zz(index) = 0;

p1  = trapz(x1,pdf1)
p2  = trapz(x2,pdf2)
p12 = trapz(x_short,zz)

plot(x1,pdf1,x2,pdf2,x12,pdf12)
hold on;
fill(x_short,zz,'blue')     % plot x1+x2
hold off;
于 2015-12-30T06:05:13.097 回答