Как рассчитать кумулятивную функцию плотности суммы случайных величин в С++?

У меня есть ряд случайных величин с известными, но, возможно, разными распределениями вероятностей (PDF). Мне нужно вычислить значение кумулятивной функции распределения (CDF) суммы этих случайных величин в определенной точке x.

Это включает в себя два шага:

1. Рассчитайте свертку PDF случайной величины, чтобы получить PDF их суммы

2. Вычислить CDF(x), т. е. интеграл суммарной PDF от -бесконечности до x

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

Например, если у меня есть одна случайная величина, распределенная по логарифмическому закону с PDF:

PDF_1(x) = e^(-(log(x) -μ)^2/(2 sigma^2))/(sqrt(2 π) x sigma)

и еще одна случайная величина, также распределенная логарифмически нормально, но с другим параметром:

PDF_2(x) = e^(-(log(x) -μ)^2/(2 delta^2))/(sqrt(2 π) x delta)

то PDF их суммы (Convolution[PDF_1, PDF_2]) не может быть рассчитан аналитически.

Предполагая произвольную точку отсечки хвоста распределения в 1000, и предполагая, что параметры μ, сигма и дельта известны:

PDF_sum(x) = NumericallyIntegrate (PDF_1(y) PDF_2(y - x) dy), y=0 to 1000

В настоящее время:

CDF(x) = NumericallyIntegrate (PDF_sum(y) dy), y=0 to x

Как написать этот код для вычисления CDF на C++, и какие библиотеки вы предлагаете?