вычисление суммы двух треугольных случайных величин (Matlab)

Я хотел бы вычислить сумму двух треугольных случайных величин,

P(x1+x2 < y)

изображение

Есть ли более быстрый способ реализовать сумму двух треугольных случайных величин в Matlab?

РЕДАКТИРОВАТЬ: кажется, что есть гораздо более простой способ, как показано в этой minitab демонстрация. Так что это не невозможно. К сожалению, это не объясняет, как был рассчитан PDF. Все еще изучаю, как я могу сделать это в Matlab.

EDIT2: Следуя совету, я использую функцию conv в Matlab для разработки 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. PDF X1 + X2 интегрируется до значения, намного превышающего 1.
  2. PDF X1+X2 широко варьируется в зависимости от диапазона x. Интуитивно понятно, что если X1+X2 больше 210 (сумма верхних пределов «с» двух отдельных треугольных распределений, 100 + 110), не должно ли P(X1+X2 ‹210) равняться 1? Кроме того, поскольку нижние пределы «а» равны 85 и 90, P(X1+X2 ‹85) = 0?

matlab statistics convolution probability probability-distribution
person Yohan K.    schedule 21.12.2015    source источник
comment
Чтобы выполнить поэлементное умножение, вам нужно .* между двумя PDF-файлами. Итак, fun = @(x) pdf(pd2,x).*pdf(pd1,y-x);   -  person schvaba986    schedule 21.12.2015
comment
Что касается терминологии, функция вероятности, которую вы определяете, зависит от одного аргумента, поэтому я бы не назвал ее «совместной» вероятностью.   -  person Itamar Katz    schedule 21.12.2015
comment
Отредактировано за неправильное понимание предмета. Извиняюсь! Теперь вопрос правильно сформулирован.   -  person Yohan K.    schedule 21.12.2015

Ответы (2)


arrow_upward
2
arrow_downward

PDF суммы независимых переменных является сверткой PDF переменных. Итак, вам нужно вычислить свертку двух переменных с помощью треугольных PDF-файлов. Треугольник кусочно-линейный, поэтому свертка будет кусочно-квадратичной.

Есть несколько способов об этом. Если численный результат приемлем: дискретизируйте PDF-файлы и вычислите свертку дискретизированных PDF-файлов. Я считаю, что для этого в Matlab есть функция conv. Если нет, вы можете выполнить быстрое преобразование Фурье (через fft), вычислить произведение по точкам, затем выполнить обратное преобразование (ifft, если я правильно помню), поскольку fft(convolution(f, g)) = fft(f) fft (грамм). Вам нужно быть осторожным с нормализацией, если вы используете conv или fft.

Если вам нужен точный результат, свертка — это просто интеграл, и если вы осторожны с пределами интегрирования, вы можете вычислить его вручную. Я не знаю, доступен ли вам символьный набор инструментов Matlab, и если да, то я не знаю, может ли он обрабатывать интегралы функций, определенных кусочно.

person Robert Dodier    schedule 21.12.2015
comment
Спасибо за ваше руководство. Кажется, conv это правильный способ сделать это в Matlab. Я обновил свой вопрос, чтобы показать прогресс. - person Yohan K.; 22.12.2015
comment
Похоже, у вас есть большая часть решения. Вы должны быть осторожны с нормализацией PDF, чтобы убедиться, что PDF x, y и x + y все интегрируются в 1. Также обратите внимание, что поддержка (ненулевая плотность) x + y - это интервал (min(x) + min(y), max(x) + max(y)), которые вы должны сопоставить с результатом conv. В частности, первый элемент результата conv соответствует первому элементу x + первому элементу y, точно так же последний элемент соответствует последнему элементу x + последнему элементу y. В результате conv есть элементы length(x) + length(y) - 1. - person Robert Dodier; 22.12.2015

arrow_upward
1
arrow_downward

Ниже приведена правильная реализация для будущих пользователей. Большое спасибо Роберту Додье за ​​руководство.

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;
person Yohan K.    schedule 30.12.2015