Задача на интегрирование из теории вероятности

Дана неотрицательная функция p(x), непрерывная на отрезке [a, b] и принимающая нулевые значения вне этого отрезка.

Дана ещё одна неотрицательная функция q(y), непрерывная на отрезке [c, d] и принимающая нулевые значения вне этого отрезка.

Мы полагаем далее, что X – непрерывная случайная величина с плотностью вероятности P ∙ p(x). Аналогично, Y – непрерывная случайная величина с плотностью вероятности Q ∙ q(y).

Требуется с помощью численного метода интегрирования высокого порядка точности вычислить с шагом по аргументу h = 0,01:

1) Постоянные множители X и Y, при которых ∫_a^b P ∙ p(x) dx = 1 и ∫_c^d Q ∙ q(y) dy = 1

2) Множество элементарных исходов случайной величины Z = X + Y

3) Плотность вероятности суммы случайных величин Z = X + Y по формуле: f(z) = ∫_-_∞^+∞ P ∙ p(x) ∙ Q ∙ q(z - x) dx, z ∈ [a + c, b + d]; 0 иначе. Пределы интегрирования станут конечными, но их следует аккуратно подобрать.

4) Также вывести график плотностей вероятности для величин X, Y, Z.

5) Полную вероятность для случайной величины Z: ∫_Zmin^Zmax f(z) dz, сравнив в конце результат с единицей.

В качестве исходных данных зададим функции: p(x) = x на отрезке [a, b] = [1, 3]; q(y) = y² на отрезке [c, d] = [0, 1].

Код:

import matplotlib.pyplot as plt
import numpy as np

#метод Ньютона для функций
def nyton(y, a, b, h):
    s = 0
    while round(a, 8) < b:
        s += (h / 8) * (y(a) + 3 * y(a + h / 3) + 3 * y(a + h / 1.5) + y(a + h))
        a += h
    return s

#метод трапеций для массива точек
def trapeze(X, Y):
    h = abs(X[1] - X[0])
    s = 0.0
    for i in range(X.shape[0] - 1):
        s += 0.5 * h * (Y[i] + Y[i + 1])
    return s

#функция p(x) на отрезке
a, b = 1, 3
def p(x):
    if (a <= x <= b):
        return x
    else:
        return 0

#функция q(y) на отрезке
c, d = 0, 1
def q(y):
    if (c <= y <= d):
        return y ** 2
    else:
        return 0

h = 0.01 #шаг

P = 1 / nyton(p, a, b, h)
Q = 1 / nyton(q, c, d, h)
print(P)
print(Q)

#множество элементарных исходов случайной величины Z=X+Y (в конкретных точках)
print(P * p(a) * Q * q(d))

#плотность вероятности суммы случайных величин Z=X+Y
n1 = (b - a) / h
n2 = (d - c) / h
#массив точек (количество n1 + n2 + 1), полученный разбиением отрезка [а + с, b + d]
Z = np.linspace(a + c, b + d, int(n1 + n2 + 1))
F = np.zeros(int(n1 + n2 + 1))
for i in range(F.shape[0]):
    F[i] = nyton(lambda x: P * p(x) * Q * q(Z[i] - x), a, b, h)

#графики плотностей вероятности
xi_p = np.linspace(a - 1, b + 1, int(n1 + 1))
yi_p = [p(i) for i in xi_p]
plt.plot(xi_p, yi_p)
plt.xlabel('x')
plt.ylabel('p(x)')
plt.show()

xi_q = np.linspace(c - 1, d + 1, int(n2 + 1))
yi_q = [q(i) for i in xi_q]
plt.plot(xi_q, yi_q)
plt.xlabel('y')
plt.ylabel('q(y)')
plt.show()

plt.plot(Z, F)
plt.xlabel('z')
plt.ylabel('f(z)')
plt.show()

#полная вероятность для случайной величины Z
print(trapeze(Z, F))

Результаты и выводы:

1) P = 0.2500000000000006, Q = 3.0112923462986156

2) Множество элементарных исходов случайной величины Z=X+Y для p(a) и q(d): 0.7528230865746558

3) -

4) Графики плотностей вероятности для величин X, Y, Z:

5) Полная вероятность для случайной величины Z: ∫_Zmin^Zmax f(z) dz = 1.002947741530758, с учетом погрешности = 1.

Статьи

Задача на интегрирование из теории вероятности