Свертка двух зависимых распределений в MATLAB - PullRequest
Новая
       17

Свертка двух зависимых распределений в MATLAB

1 голос
/ 15 мая 2019

Предположим, что у меня есть две дискретные случайные величины X и Y . 

X = {1,3,3,5,7,7,7,9,9,9,9,9}

и 

Y = {5,5,9,9,10,12,13}

Где их эмпирические CDF даны как: 

F_x(1) = 0.0833, F_x(3) = 0.25, F_x(5) = 0.33, F_x(7) = 0.5833 and F_x(9) = 1

и 

F_y(5) = 0.2857, F_y(9) = 0.5714, F_y(10) = 0.7143, F_y(12) = 0.8571 and F_y(13) = 1

Предполагая, что их совместное распределение 

H(x,y) = F_x(x) * F_y(y)

что на самом деле является "допущением" X и Y независимы.

Как можно вычислить Z = X + Y и F (z) в MATLAB?

Примечание: я дал H (x, y) как простую функцию произведения для простоты, но в действительности это может быть что угодно, что фактически моделирует зависимость между X и Y.

1 Ответ

1 голос
/ 15 мая 2019

С учетом непрерывных функций плотности вероятности F X и F Y с совместной функцией плотности вероятности F X, Y , мы можем вычислить F X + Y как интеграл от F X, Y по линии z = x + y . Если функции плотности вероятности являются дискретными, вышеприведенный интеграл должен быть записан как производная интеграла по части плоскости, заданная как z <= <em>x + y .

Это довольно просто сделать в MATLAB. Начнем с данных ОП: 

F_x = [0.0833,0.0833,0.25,0.25,0.33,0.33,0.5833,0.5833,1]; % CDF
F_x = diff([0,F_x]); % PDF
F_y = [0,0,0,0,0.2857,0.2857,0.2857,0.2857,0.5714,0.7143,0.7143,0.8571,1]; % CDF
F_y = diff([0,F_y]); % PDF
H = F_x.' .* F_y; % example joint PDF

Теперь мы суммируем F_cum(z) = sum(H(x,y)) для всех значений z<=x+y, а затем берем производную F = diff([0,F_cum]): 

[m,n] = size(H);
F_cum = zeros(1,m+n-1);
for z = 1:numel(F_cum)
   s = 0;
   for x = 1:numel(F_x)
      y = z-x+1;
      y = max(min(y,n),1); % avoid out of bounds indexing
      s = s + sum(H(x,1:y));
   end
   F_cum(z) = s;
end
F = diff([0,F_cum]);

Обратите внимание, что мы определили y=z-x+1, что означает z=y+x-1. Таким образом, F(1) соответствует z=2. Это наименьшее возможное значение, которое может получиться из суммы двух распределений, которые мы определили, чтобы начать с 1.

Вышесказанное можно упростить, добавив к нулю H и сдвинув каждую строку на один дополнительный элемент. Это выстраивает линию z = x + y в столбце матрицы, что позволяет нам использовать тривиальную проекцию sum: 

H = [H,zeros(m)];
for ii=2:m
   H(ii,:) = circshift(H(ii,:),ii-1);
end
F_cum = cumsum(sum(H,1));
F_cum = F_cum(1:end-1); % last element we don't need
F2 = diff([0,F_cum]);

Но поскольку diff([0,cumsum(F)]) == F (с точностью до чисел), мы можем пропустить эти две операции: 

F3 = sum(H,1);
F3 = F3(1:end-1); % last element we don't need

(all(abs(F-F2)<1e-15) и all(abs(F-F3)<1e-16))

...