1

2 つの三角確率変数の合計を計算したいのですが、

P(x1+x2 < y)

画像

Matlab で 2 つの三角確率変数の合計を実装するより高速な方法はありますか?

編集: このminitabデモンストレーションに示されているように、もっと簡単な方法があるようです。だから、それは不可能ではありません。残念ながら、PDFがどのように計算されたかは説明されていません。matlabでこれを行う方法をまだ調べています。

EDIT2: アドバイスに従って、convMatlab の関数を使用して、2 つの確率変数の合計の 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;

ただし、このコードには 2 つの問題があります。

  1. X1+X2 の PDF を統合すると、1 よりもはるかに高くなります。
  2. X1+X2 の PDF は、x の範囲によって大きく異なります。直観的に、X1+X2 が 210 (2 つの個別の三角分布の上限 "c" の合計、100 + 110) より大きい場合、P(X1+X2 <210) は 1 に等しくないでしょうか? また、下限ａは８５と９０なので、Ｐ（Ｘ１＋Ｘ２＜８５）＝０？
4

2 に答える 2

2

独立変数の合計の pdf は、変数の pdf の畳み込みです。したがって、三角pdfを使用して2つの変数の畳み込みを計算する必要があります。三角形は区分線形であるため、畳み込みは区分二次になります。

それについてはいくつかの方法があります。数値結果が許容できる場合: pdf を離散化し、離散化された pdf の畳み込みを計算します。そのための関数convがMatlabにあると思います。そうでない場合は、（を介して）高速フーリエ変換を実行し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 に答える