Метод кусочной аппроксимации плотности распределения вероятности

ГУАП

КАФЕДРА № 21

ОТЧЕТ
ЗАЩИЩЕН С ОЦЕНКОЙ

ПРЕПОДАВАТЕЛЬ

РАБОТУ ВЫПОЛНИЛ

Санкт-Петербург 2025

1. Цель

В соответствии с номером варианта необходимо сгенерировать выборку случайных чисел с заданным законом распределения с помощью метода Брусленко Н.П., исследовать данный метод.

1. Выполнение работы

 Заданный закон распределения[pic 1]

1. Подставляем известные значения[pic 2][pic 3]

Рисунок 1 – График, построенный по формуле плотности распределения

1. Алгоритм

Алгоритм работы программы поделен мной на 6 блоков:

      - для подпрограмм:

А)   создание функции с заданным законом распределения;

Б) создание функции численного интегрирования методом средних прямоугольников;

- для основной программы:

В) Численный поиск максимума функции;

Г) Поиск правой границы b (где функция < 1% от максимума);

Д) Подготовка данных для метода Брусленко;

Е) Генерация случайных чисел методом Брусленко[pic 4]

Рисунок 1 – Алгоритм программы блока А

                       а)                                                                           б)[pic 5][pic 6]

Рисунок 2 (а,б) – Алгоритм работы программы: а) блока Б; б) блока В[pic 7]

 Рисунок 3 – Алгоритм работы программы блока Г

                             а)                                                                         б) [pic 8][pic 9]

Рисунок 4 (а,б) – Алгоритм работы программы: а) блока Д, б) блока Е

1. Программа (C++)

1. #include <iostream>
2. #include <cmath>
3. #include <ctime>
4. #include <cstdlib>
6. using namespace std;
7. /*
8. Функция плотности вероятности
9. y - входное значение
10. Возвращает значение функции плотности в точке y
11. */
12. double func(double y) {
13.     return  pow(y, 3) * exp(-y * y / 2.0) / 2.0;
14. }
15. /*
16.  * Функция численного интегрирования методом средних прямоугольников
17.  * a - нижний предел интегрирования
18.  * b - верхний предел интегрирования
19.  * Возвращает значение интеграла от a до b с точностью E
20.  */
21. double integ(double a, double b) {
22.     const double E = 0.001; // желаемая точность вычисления интеграла
23.     int n = 1; // начальное количество интервалов разбиения
24.     double h = (b - a) / n; // Шаг интегрирования
25.     double S1 = 0, S2; // S1 - текущая сумма, S2 - предыдущая сумма

1.     do {
2.         n *= 2; // Удваиваем количество интервалов
3.         h = (b - a) / n; // Пересчитываем шаг интегрирования
4.         S2 = S1; // Сохраняем предыдущее значение суммы
5.         S1 = 0; // Обнуляем текущую сумму
7.         // Суммируем значения функции в средних точках интервалов
8.         for (int i = 0; i < n; ++i) {
9.             S1 += func(a + (i + 0.5) * h);
10.         }
12.         S1 *= h; // Умножаем сумму на шаг интегрирования
14.         // Повторяем, пока не достигнем нужной точности
15.     } while (fabs(S1 - S2) > E);
16.     return S1;
17. }
19. int main() {
20.     // Инициализация генератора случайных чисел текущим временем
21.     srand(static_cast<unsigned>(time(nullptr)));
23.     // Численный поиск максимума функции
24.     double max_y = 0;       // Максимальное значение функции
25.     double x = 0;       // Точка, где достигается максимум
26.     double s = 0.001;    // Шаг поиска
27.     double y0 = 0.1; // Начальная точка поиска (избегаем y=0)
29.     // Ищем максимум на интервале [0.1, 10.0]
30.     while (y0 <= 10.0) {
31.         double val = func(y0); // Значение функции в текущей точке
33.         // Если нашли большее значение, обновляем максимум
34.         if (val > max_y) {
35.             max_y = val;
36.             x = y0;
37.         }
39.         // Переходим к следующей точке
40.       y0 += s;
41.     }
42.     cout << "max y = " << max_y << endl;
44.     // Поиск правой границы b
45.     double d = 0.1; // Шаг поиска границы
46.     double y = x; // Начинаем поиск от точки максимума
47.     double b = x; // Правая граница области
49.     // Двигаемся вправо, пока значение функции > 1% от максимума
50.     while (func(y) > max_y * 0.01) {
51.         y += d;
52.         b = y;
53.     }
54.     cout << "b = " << b << endl;
56.     // Подготовка данных для метода Брусленко
57.     const int K = 100; // Количество интервалов разбиения
58.     double yk[K + 1];  // Массив границ интервалов
59.     double F[K + 1] = { 0 };  // Массив значений интегральной функции
61.     // Равномерное разбиение интервала [0, b] на K частей
62.     for (int k = 0; k <= K; ++k) {
63.         yk[k] = k * b / K; // Вычисляем k-ю границу
64.     }
66.     // Вычисление интегральной функции распределения
67.     for (int k = 1; k <= K; ++k) {
68.         F[k] = F[k - 1] + integ(yk[k - 1], yk[k]);
69.     }
71.     //Нормировка интегральной функции
72.     for (int k = 1; k <= K; ++k) {
73.         F[k] /= F[K]; // Делим все значения на F[K] (интеграл по всей области)
74.     }
76.     // Генерация случайных чисел
77.     for (int i = 0; i < 100; ++i) {
78.         //Генерируем равномерно распределенное число c1 ∈[0, 1]
79.         double c1 = static_cast<double>(rand()) / RAND_MAX;
80.         //Находим интервал k, где F[k] <= c1 < F[k + 1]
81.         int k = 0;
82.         while (k < K && F[k + 1] < c1) {
83.             k++;
84.         }
85.         //Генерируем второе равномерное число c2 ∈ [0,1]
86.         double c2 = static_cast<double>(rand()) / RAND_MAX;
88.         //Линейная интерполяция внутри выбранного интервала
89.         double n = yk[k] + (yk[k + 1] - yk[k]) * c2;
90.         cout << n << endl;
91.     }
93.     return 0;
94. }

                                   а)                                                                                                   б)[pic 10][pic 11]