计算两个三角随机变量的总和(Matlab)
calculating sum of two triangular random variables (Matlab)
我想计算两个三角随机变量的总和,
P(x1+x2 < y)
有没有更快的方法在Matlab中实现两个三角随机变量之和?
编辑:看起来可能有更简单的方法,如 minitab 演示所示。所以这不是不可能的。遗憾的是,它没有解释 PDF 是如何计算的。仍在研究如何在 matlab 中执行此操作。
EDIT2:根据建议,我在 Matlab 中使用 conv
函数来开发两个随机变量之和的 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;
但是这段代码有两个问题:
- X1+X2 的 PDF 集成度远高于 1。
- X1+X2 的 PDF 根据 x 的范围变化很大。直观上,如果X1+X2大于210(两个个体三角分布的上限"c"之和,100+110),P(X1+X2<210)不应该等于1吗?另外,由于下限 "a" 是 85 和 90,P(X1+X2 <85) = 0?
自变量之和的 pdf 是变量 pdf 的卷积。所以你需要用三角 pdf 计算两个变量的卷积。三角形是分段线性的,因此卷积将是分段二次的。
有几种方法可以解决这个问题。如果数值结果是可接受的:将 pdf 离散化并计算离散化 pdf 的卷积。我相信在 Matlab 中有一个函数 conv
。如果没有,您可以采用快速傅里叶变换(通过 fft
),逐点计算乘积,然后进行逆变换(ifft
,如果我没记错的话)因为 fft(convolution(f, g )) = fft(f) fft(g)。如果您使用 conv
或 fft
.
,则需要注意规范化
如果你一定要有一个精确的结果,卷积只是一个积分,如果你注意积分的极限,你可以手算出来。不知道Matlab的symbolic toolbox是否对你可用,如果有,我不知道它是否可以处理分段定义的函数的积分。
以下是适合未来用户的正确实施。非常感谢 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;
我想计算两个三角随机变量的总和,
P(x1+x2 < y)
有没有更快的方法在Matlab中实现两个三角随机变量之和?
编辑:看起来可能有更简单的方法,如 minitab 演示所示。所以这不是不可能的。遗憾的是,它没有解释 PDF 是如何计算的。仍在研究如何在 matlab 中执行此操作。
EDIT2:根据建议,我在 Matlab 中使用 conv
函数来开发两个随机变量之和的 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;
但是这段代码有两个问题:
- X1+X2 的 PDF 集成度远高于 1。
- X1+X2 的 PDF 根据 x 的范围变化很大。直观上,如果X1+X2大于210(两个个体三角分布的上限"c"之和,100+110),P(X1+X2<210)不应该等于1吗?另外,由于下限 "a" 是 85 和 90,P(X1+X2 <85) = 0?
自变量之和的 pdf 是变量 pdf 的卷积。所以你需要用三角 pdf 计算两个变量的卷积。三角形是分段线性的,因此卷积将是分段二次的。
有几种方法可以解决这个问题。如果数值结果是可接受的:将 pdf 离散化并计算离散化 pdf 的卷积。我相信在 Matlab 中有一个函数 conv
。如果没有,您可以采用快速傅里叶变换(通过 fft
),逐点计算乘积,然后进行逆变换(ifft
,如果我没记错的话)因为 fft(convolution(f, g )) = fft(f) fft(g)。如果您使用 conv
或 fft
.
如果你一定要有一个精确的结果,卷积只是一个积分,如果你注意积分的极限,你可以手算出来。不知道Matlab的symbolic toolbox是否对你可用,如果有,我不知道它是否可以处理分段定义的函数的积分。
以下是适合未来用户的正确实施。非常感谢 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;