Зона кода

А давайте немного попрограммируем...

Генерация нормально распределённых псевдослучайных чисел

C99Решение задач

Просматривая книгу Седжвика, Уэйна и Дондеро "Программирование на языке Python: учебный курс" (я уже упоминал её на страницах данного сайта здесь и здесь), я наткнулся на код программы "Броуновский мост". В ней вызывалась функция gaussian(), входящая в стандартный модуль stdrandom языка Python.

Эта функция каждый раз при её вызове возвращает число из генерируемой ею последовательности псевдослучайных чисел, распределённых по нормальному закону. И тут я задумался о том, как создать аналог данной функции на языке C99.

Ведь все стандартные функции языка C99, каким-либо образом связанные с генерацией псевдослучайных чисел, исчерпываются всего лишь двумя функциями, входящими в библиотеку stdlib.h. Одна из них — это rand(). Данная функция генерирует псевдослучайные целые числа, равномерно распределённые в диапазоне от 0 до RAND_MAX включительно (в нашем случае значение макроса RAND_MAX равно 32767).

Другая функция — это srand(), которая "вынуждает" функцию rand() генерировать последовательность, "жёстко связанную" с аргументом, передаваемым srand() (каждому такому аргументу соответствует своя последовательность).

Вот и все функции, которыми мы располагаем. Я в данной статье решаю задачу построения функции на языке C99, генерирующей последовательность псевдослучайных чисел, распределённых по нормальному закону. Для понимания статьи читателю пригодятся знания теории вероятностей в том объёме, в котором она даётся в техническом вузе. Также нелишним будет владение элементарными методами вычислительной математики.

Выбор варианта решения поставленной задачи

Напомню, что случайная величина, распределённая по нормальному закону, — это непрерывная случайная величина, заданная следующей функцией распределения:

Fm,σx=1σ2π-xe-t-m22σ2dt.

Здесь m и σ — параметры случайной величины, имеющие смысл её математического ожидания и среднеквадратического отклонения соответственно. Часто нормальное распределение называют также распределением Гаусса.

В дальнейшем мы будем рассматривать стандартную нормальную случайную величину, имеющую параметры m = 0 и σ = 1. Её функцию распределения обозначим F(x):

Fx=12π-xe-t22dt.

Стоит нам научиться генерировать псевдослучайные числа, распределённые по стандартному нормальному закону, так сразу же мы сможем получить и числа, распределённые по нормальному закону с произвольными значениями параметров (об этом поговорим в конце статьи).

Заметим, что с помощью функции rand() можно сгенерировать псевдослучайные числа, равномерно распределённые на отрезке [0; 1]. Для этого достаточно каждое значение, возвращаемое rand(), делить на значение макроса RAND_MAX. Так что постановку нашей задачи можно конкретизировать и перевести в теоретическую плоскость, не учитывающую псевдослучайность значений, получаемых программно.

В этой постановке задача будет формулироваться так: получить случайную величину, распределённую по стандартному нормальному закону, имея в распоряжении непрерывную случайную величину, равномерно распределённую на отрезке [0; 1]. Сразу же замечу, что эта задача имеет достаточно простое решение. Обозначим, например, первую из упомянутых величин символом Y, а вторую — символом X. Легко показать, что Y можно выразить через X следующим образом:

Y = F −1(X).

Здесь F −1(x) — это функция, обратная по отношению к функции F(x) (такая функция существует, поскольку F(x) монотонна на всей числовой оси).

Основная проблема при практическом использовании такого решения заключается в сложности вычисления значений функции F −1(x).

Имеются и другие способы получения нормально распределённой случайной величины из распределённой равномерно. Некоторые из них описаны во 2-м томе "Искусства программирования" Дональда Кнута (способ, приведённый выше, там же тоже рассмотрен).

Из алгоритмов, приведённых в книге Кнута, самыми простыми для реализации являются, пожалуй, метод полярных координат и метод отношений. Оба этих метода используют в качестве исходных две независимые непрерывные случайные величины, распределённые равномерно на отрезке [0; 1], причём первый из них приводит к получению "на выходе" сразу двух нормально распределённых независимых случайных величин.

Ознакомившись с алгоритмами, предлагаемыми Кнутом, я решил, всё же, остановиться на самом первом варианте, приведённом в данном разделе, сводящем решение задачи к вычислению F −1(x). Объясню причины такого решения.

Во-первых, даже упомянутые мною два алгоритма из книги Кнута требуют наличие "на входе" двух независимых равномерно распределённых на отрезке [0; 1] случайных величин. Если мы сразу же задумаемся о программной реализации этих алгоритмов, то перед нами сразу встанет вопрос: а сможем ли мы смоделировать две такие независимые случайные величины, обладая лишь одним единственным генератором псевдослучайных чисел (функцией rand())?

Ведь в нашем распоряжении имеется только одна последовательность псевдослучайных чисел. Допустим, мы разделим её на 2 подпоследовательности (по принципу чётности и нечётности членов). В какой степени эти две новые последовательности можно рассматривать как модели независимых случайных величин, если они получены из одной последовательности? На этот вопрос я ответить не готов; для ответа, как мне кажется, требуется проводить дополнительные исследования.

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

Вторая причина, побудившая меня обратиться к методу, требующему вычисления F −1(x), заключается в том, что стандартная библиотека языка C99 может мне "помочь" с вычислением функции F (x). К сожалению, вычисление значений самой функции F (x) в библиотеке не реализовано, но зато библиотекой math.h поддерживается вычисление функции erf(x), называемой функцией ошибок, выражающейся через интеграл следующим образом:

erfx=2π0xe-t2dt.

Вычислением её значений занимается функция erf() из библиотеки math.h. Подчеркну, что поддержка данной функции предусмотрена стандартом C99, но стандартные библиотеки более ранних версий языка C поддерживать эту функцию не обязаны.

Так вот, функция erf(x) является "родственной" по отношению к F(x). Говоря более конкретно: F(x) можно выразить через erf(x) следующим образом:

F(x)=121+erfx2.

А для чего нам нужно уметь вычислять значения F(x), ведь нас же интересует функция F−1(x)? Ответ на этот вопрос читатель получит, ознакомившись со следующим разделом статьи.

Вычисление значений функции F −1(x)

Предположим, что мы хотим вычислить значение функции F −1(x) в некоторой точке x = a (разумеется, число a предполагаем принадлежащим множеству определения функции F −1(x), т. е. интервалу (0; 1)). Искомое значение F −1(a), очевидно, является корнем уравнения

F(x) = a.

Ниже приведён рисунок, иллюстрирующий данный факт. Он содержит графики функций y = F(x) и y = a. Видно, что абсцисса точки пресечения графиков — это и есть F −1(a).

Графики функций y = F(x) и y = a

Графики функций y = F(x) и y = a (щёлкните для увеличения)

Давайте составим схему приближённого решения уравнения F(x) = a одним из численных методов. Я предлагаю применить весьма простой в реализации метод, называемый методом Ньютона или методом касательных. Я уже упоминал и даже описывал этот метод в статьях данного сайта (здесь и здесь), так что на этот раз буду краток.

Метод Ньютона позволяет приближённо решать уравнения вида

g(x) = 0,

при условии, что функция g(x) имеет в некоторой окрестности X корня x~ этого уравнения отличную от нуля производную. Строится задаваемая рекуррентно последовательность следующего вида:

x0 задано, xk+1=xk-gxkg'xk, k=0,1,2, .

Здесь x0 — это некоторое начальное значение, принадлежащее X. Разумеется, встаёт вопрос о существовании членов этой последовательности. Однако, известно, что если все они существуют и принадлежат X, то данная последовательность сходится к корню исходного уравнения.

На практике, естественно, члены последовательности прекращают вычислять при достижении нужной степени точности и последний вычисленный член берут в качестве приближённого корня уравнения.

В качестве условия окончания итераций (т. е. вычислений членов последовательности) часто выбирают выполнение неравенства

|xk+1xk| < ε,

где ε — некоторое заранее выбранное малое положительное число. При этом рассчитывают на то, что абсолютная погрешность приближённого значения xk+1 (т. е. x~-xk+1) не превысит ε. Такой же критерий окончания итераций будем использовать и мы.

В нашем случае, очевидно, g(x) = F(x) − a. Таким образом, итерационная схема будет иметь вид.

x0 задано, xk+1=xk-Fxk-aF'xk, k=0,1,2, .

Учитывая связь между F(x) и erf(x), а также то, что

F'x=12πe-x22,

перепишем схему в виде

x0 задано, b=1-2a, xk+1=xk-π2erfx2+bex22, k=0,1,2, .

Как известно, функция F(t) обладает свойством

F(t) + F(−t) = 1,

из которого следует, что

F(−t) = 1 − F(t) .

Возьмём от обеих частей этого равенства функцию F −1:

t = F −1(1 − F(t)),

откуда

t = −F −1(1 − F(t)).

Сделаем замену: x = F(t), т. е. t = F −1(x). В итоге, получаем:

F −1(x) = −F −1(1 − x).

Таким образом, мы можем ограничиться реализацией вычислений значений функции F −1(x) лишь на интервале (0,5; 1), а её значения на интервале (0; 0,5) получать по приведённой выше формуле (а то, что F(0,5) = 0 мы знаем и так, без всяких вычислений).

От чего зависит количество итераций, обеспечивающих нахождение приближённого значения корня с заданной степенью точности? От двух вещей.

Во-первых, от выбора начального значения x0. Разумеется, чем ближе x0 к корню уравнения, тем меньшее количество итераций потребуется.

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

К сожалению, график функции F(x) при увеличении x "очень быстро" начинает "прижиматься снизу" к прямой y = 1 (см. рисунок), т. е. производная F(x) (совпадающая с производной функции F(x) − a) становится слабо отличимой от нуля. Это неблагоприятным образом сказывается на количестве итераций при значениях a, близких к 1. Можно попытаться компенсировать увеличение итераций при таких a выбором начальных значений x0.

Поговорим, как раз, о выборе значений x0. Поступим следующим образом. Найдём приближённые значения функции F −1(x) в "эталонных" точках 0,5, 0,6, 0,7, 0,8, 0,9, 0,97. Для приближённого нахождения значения F −1(x) в точке a в качестве начального значения x0 возьмём значение F −1(x) в той из этих шести точек, к которой точка a наиболее близка (исключение составит лишь случай попадания a в интервал (0,935; 0,95); здесь в качестве x0 будет выбрано ≈F −1(0,95)).

А что, если нам потребуется найти F −1(0) и F −1(1)? С формальной точки зрения, F −1(x) в точках 0 и 1 не определена, однако имеет в этих точках следующие односторонние пределы:

limx-0F-1x=-, limx1+0F-1x=+.

Тем не менее, случайная величина, непрерывно распределённая на отрезке [0; 1], может принимать значения, совпадающие с границами отрезка (правда, теоретически вероятности принятия таких значений нулевые, но наша случайная величина, моделируемая программно, может принять эти значения пусть и с малыми вероятностями, но отличными от нуля). Так что решать проблему "вычисления" F −1(0) и F −1(1), всё же, придётся.

Можно, конечно, возвращать специальные значения типа double, соответствующие "−∞" и "+∞", но это не лучшее решение, поскольку конструируемый нами генератор псевдослучайных чисел, распределённых по нормальному закону, должен всегда выдавать именно вещественные числа, и ничего кроме них.

Итак, какое значения приписать F −1(1)? Напомню, что мы моделируем непрерывную случайную величину, равномерно распределённую на отрезке [0; 1], вызывая функцию rand() и деля возвращённое ею значение на RAND_MAX, т. е. на 32767. Максимальное значение такого частного, не превышающее единицу, равно, очевидно, 32766 / 32767, т. е. 0,99997… . Так что F −1(1) припишем значение, немного превышающее F −1(0,99997) ≈ 4,012811. Я решил принять это значение равным 4,3. В этом случае, очевидно, в качестве F −1(0) следует брать число −4,3.

Наша функция, написанная на C99, будет возвращать значение 4,3 даже в более общем случае: если переданный ей аргумент превышает разность 1 − ε. А если значение меньше ε, то функция будет возвращать −4,3.

Какое число принять за ε? Два соседних значения, выдаваемых нашим генератором псевдослучайных чисел, равномерно распределённых на отрезке [0; 1], различаются на 1 / 32767, т. е. примерно на 3,05 * 10-5. Нет смысла вычислять значения F −1(x) с точностью, сильно превышающей точность, с которой "выдаются" упомянутые выше псевдослучайные числа. Так что в качестве ε я решил взять, для простоты, 10−5.

Ну что ж, все теоретические проблемы мы обсудили, так что можно приступать к программированию.

Создание функции, вычисляющей F −1(x) и её тестирование

Мы создадим функцию, вычисляющую F −1(x), в рамках программы, которая будет эту функцию тестировать. Вся программа целиком уместится в одном-единственном файле.

В начале этого файла размещаем инструкции, подключающие стандартные библиотеки:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

Следующий блок кода содержит определения макросов-констант:

#define sqrt_2 1.414213562373       //sqrt(2)
#define sqrt_pi_2 1.253314137316    //sqrt(pi / 2)
#define eps 1e-5                    //погрешность

Смысл первых двух констант, полагаю, ясен из комментариев. Что касается макроса-константы eps, то он имеет тот же смысл, что и число ε.

Объявляем и инициализируем массив, который содержит приближённые значения функции F −1(x) в 6 точках, которые будут использоваться в качестве начальных значений для итераций:

const double val[] = {0.0,             //f_1(0.5)
                      0.253347,        //f_1(0.6)
                      0.524401,        //f_1(0.7)
                      0.841621,        //f_1(0.8)
                      1.281552,        //f_1(0.9)
                      1.880794};       //f_1(0.97)

Объявляем глобальную переменную it — счётчик итераций:

int it;

Наконец, создаём функцию f_1(), вычисляющую значения F −1(x):

 1.double f_1(double x)
 2.{
 3.    if (x == 0.5)
 4.        return 0.0;
 5.    if (x < eps)
 6.        return -4.3;
 7.    if (x > 1 - eps)
 8.        return 4.3;
 9.    int sign = 1;
10.    if (x < 0.5)
11.    {
12.        x = 1 - x;
13.        sign = -1;
14.    }
15.    double xk1 = val[(int) (round(10 *x) - 5)], xk;
16.    double b = 1 - 2 * x;
17.    it = 0;
18.    do
19.    {
20.        xk = xk1;
21.        double t = xk / sqrt_2;
22.        xk1 = xk - sqrt_pi_2 * exp(t * t) * (erf(t) + b);
23.        it++;
24.    }
25.    while (fabs(xk - xk1) >= eps);
26.    return sign == 1 ? xk1 : -xk1;
27.}

Функция f_1() вычисляет и возвращает приближённое значение функции F −1(x) в точке, переданной f_1() в виде параметра. Вычисление производится с теми оговорками, касающимися значений аргументов, близких к 0 и 1, о которых шла речь в предыдущем разделе.

Сначала обрабатываем частные случаи аргументов функции F −1(x), не требующие привлечение метода Ньютона (стр. 3-8).

Далее (см. стр. 9-14) выясняем, в какой из интервалов попадает аргумент: в (0; 0,5) или в (0,5; 1). В первом случае от аргумента x переходим к аргументу 1 - x. В переменной sign сохраняем знак искомого значения функции F −1(x) (т. е. число −1, если знак отрицателен или число 1, если положителен).

Создаём переменные xk1 и xk — аналоги текущего (вычисляемого) и предыдущего членов последовательности приближённых значений корня соответственно; первой из этих переменных сразу присваиваем соответствующее аргументу x начальное значение из массива val (стр. 15). Обнуляем счётчик итераций (стр. 17).

Сами итерации реализованы в цикле do while (стр. 18-25). Последний вычисленный член последовательности, содержащийся в xk1, присваиваем переменной xk (стр. 20), после чего вычисляем новое значение xk1 через старое, хранящееся в xk, в соответствии со схемой, приведённой ранее, т. е. выполняем итерацию (стр. 21-22). Не забываем инкрементировать счётчик итераций (стр. 23).

Цикл продолжается, пока требуемая погрешность не превышает модуля разности двух последних вычисленных членов последовательности (см. стр 25). Как только это условие нарушается, т. е., как только модуль разности становится меньше погрешности, итерации заканчиваются; последний вычисленный член последовательности в этот момент содержится в переменной xk1.

Остаётся выполнить последнее действие: в качестве окончательного результата вернуть либо значение переменной xk1, либо значение, противоположное данному, в зависимости от знака результата, хранящегося в переменной sign (стр. 26).

Разумеется, окончательная версия функции f_1() не должна содержать кода, работающего со счётчиком итераций, т. е. строк 17 и 23. Данный код понадобится нам только для тестирования.

Собственно, количество итераций, требующихся для вычислений значений F −1(x) в различных точках, — это, как раз, то, что будет нас интересовать в ходе тестирования в первую очередь.

Создадим, для удобства, функцию calc(), выводящую на печать переданную ей в качестве параметра точку, значение функции F −1(x) в этой точке, полученное, естественно, с помощью вызова f_1(), и количество итераций, потребовавшееся для вычисления этого значения:

void calc(double x)
{
    printf("%g -> %.6f ", x, f_1(x));
    printf("(итераций: %d)\n", it);
}

Не могу упустить возможность задать вопрос новичкам в программировании (если вдруг они читают эту статью). Почему нельзя было ограничиться единственным вызовом функции printf()? Другими словами, почему тело функции calc() не может состоять из единственной строчки, приведённой ниже?

    printf("%g -> %.6f (итераций: %d)\n", x, f_1(x), it);

Ну что, догадались?

А вот и ответ.

Всё дело в способе передачи фактических параметров функции в момент её вызова. Они передаются через стек. При этом, стандарт языка не оговаривает, в каком порядке параметры в него заталкиваются.

Предположим теперь, что параметры заталкиваются в стек в обратном порядке (т. е. справа налево). Разумеется, вычисляются параметры непосредственно перед заталкиванием (нет смысла их сначала вычислять и где-то временно хранить, а затем заталкивать). Таким образом, при вызове printf() (см. приведённую выше строку кода) сначала в стек будет помещено значение переменной it, затем будет произведён вызов f_1(x)и в стек будет направлен его результат, далее в стек будет помещено значение x, и, наконец, последним в стек пойдёт адрес C-строки, описывающей формат вывода.

Таким образом, в стеке окажется (и затем будет выведено на печать) то значение переменной it, которое ещё не изменено вызовом функции f_1(), т. е. "старое" значение. А это совсем не то, что нам нужно!

Что касается двух вызовов функции printf() (см. код функции calc()), то они гарантируют, что значение it будет выведено на печать уже после вызова f_1().

Вообще, программистам рекомендуется не делать никаких предположений о способе помещения в стек фактических параметров функции, т. е. писать универсальный код, корректный при любом способе. Кстати, я, в своё время, проводил эксперимент, который показал, что в программах, построенных с помощью компилятора MinGW64, параметры заталкиваются в стек именно в обратном порядке. Об этом эксперименте можно почитать здесь.

И ещё. Проблема обсуждаемого вызова функции printf() обусловлена ещё и тем, что функция f_1() обладает так называемым побочным эффектом: она изменяет значение глобальной переменной it. Такого рода побочные эффекты часто ведут к ошибкам, которые не так уж и прост обнаружить. Именно по этой причине часто рекомендуют свести использование глобальных переменных к минимуму (а, по возможности, не использовать их вовсе).

Но нам, для подсчёта количества итераций удобно счётчик сделать, всё-таки, глобальным.

Очень уж долгим оказалось отступление! Возвращаемся к нашей тестирующей программе. Непосредственно тестированием будет заниматься функция main(), имеющая следующий незатейливый код:

int main()
{
    calc(0.55);
    calc(0.65);
    calc(0.75);
    calc(0.85);
    calc(0.94);
    calc(0.95);
    calc(0.98);
    calc(0.99);
    calc(0.999);
    calc(0.9999);
    calc(0.99997);
    calc(0.99999);
    return 0;
}

Мы просто берём и вычисляем приближённые значения F −1(x) в различных точках и выясняем, сколько требуется итераций для каждого такого вычисления. Несколько первых точек специально взяты как средние арифметические между "эталонными" точками, приближённые значения F −1(x) в которых хранятся в массиве val. Таким образом мы пытаемся в некоторой степени минимизировать влияние близости выбранных нами точек к "эталонным" на количества итераций.

Вот результат выполнения нашей программы:

0.55 -> 0.125661 (итераций: 3)
0.65 -> 0.385320 (итераций: 3)
0.75 -> 0.674490 (итераций: 4)
0.85 -> 1.036433 (итераций: 4)
0.94 -> 1.554774 (итераций: 4)
0.95 -> 1.644854 (итераций: 4)
0.98 -> 2.053749 (итераций: 4)
0.99 -> 2.326348 (итераций: 5)
0.999 -> 3.090232 (итераций: 8)
0.9999 -> 3.719016 (итераций: 10)
0.99997 -> 4.012811 (итераций: 11)
0.99999 -> 4.264891 (итераций: 12)

Итак, что мы видим? Как и ожидалось, при увеличении аргумента функции F −1(x), количество итераций, необходимых для нахождения её приближённого значения, возрастает. Действительно, ведь увеличение аргумента ведёт к приближению к 0 производной функции F(x) − a, чего очень "не любит" метод Ньютона.

4 итерации — это вполне нормально, на мой взгляд. А вот 12 — это уже очень много (впрочем, как мы помним, при использовании нашего компилятора, делением значения, возвращаемого rand(), отличного от RAND_MAX, на RAND_MAX, мы можем получить максимум 0,99997… (а данному числу соответствуют лишь 11 итераций), но никак не 0,99999)! Как мы видим, выход за предел 4-х итераций наступает после превышения аргументом числа 0,98.

Можно, конечно же, попытаться число итераций уменьшить, увеличив количество "эталонных" точек в самом конце интервала (0; 1). Но стоит ли усложнять функцию ради этого? Я считаю, что нет. Ведь (псевдо)случайная величина, равномерно распределённая на отрезке [0; 1], принимает значения, превышающие 0,98 или меньшие 0,02 достаточно редко. А при генерации большого количества псевдослучайных чисел, распределённых по нормальному закону, важную роль играет среднее число итераций, а не максимальное.

И если среднее число итераций не превышает четырёх, то и улучшать нашу функцию, думаю, не стоит. А подсчётом среднего числа итераций мы ещё будем заниматься.

Создание функции, генерирующей псевдослучайные числа, распределённые по нормальному закону, и её тестирование

Функция get_norm_rand_val() будет генерировать псевдослучайные числа, распределённые по стандартному нормальному закону. По сути, она будет преобразовывать одну последовательность псевдослучайных чисел, генерируемую функцией rand(), в другую. Вот код этой функции:

double get_norm_rand_val()
{
    return f_1(rand() / ((double) RAND_MAX));
}

Думаю, смысл кода настолько очевиден, что не требует комментариев.

Нужно иметь в виду, что, многократно вызывая функцию get_norm_rand_val() (а каким образом её ещё можно использовать?), следует не вызывать в промежутках между её вызовами функцию rand(). Такие промежуточные вызовы могут привести к тому, что числа, полученные вызовами rand() внутри get_norm_rand_val(), уже не будут распределены равномерно на отрезке [0; 1], и, как следствие, числа, генерируемые get_norm_rand_val(), не будут распределены по нормальному закону.

Итак, каким образом мы будем производить тестирование нашего генератора нормально распределённых псевдослучайных чисел? Разумеется, сначала мы сгенерируем некий набор таких чисел — выборку. Более или менее серьёзное исследование такой выборки подразумевает проверку гипотезы о нормальном распределении чисел, образующих выборку, например, с помощью критерия согласия Пирсона (χ2-критерия).

Но мы пойдём гораздо более простым путём: построим гистограмму относительных частот и наложим на неё график плотности распределения вероятностей. Затем качественно, "на глаз", оценим, в какой степени гистограмма "соответствует" графику.

Напомню, что плотность распределения вероятностей, — это, ни что иное, как F '(x).

А ещё не забудем подсчитать количество всех итераций, использовавшихся для формирования выборки, чтобы затем вычислить среднее число итераций, приходящееся на однократный вызов функции get_norm_rand_val().

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

Итак, приступаем!

Модернизируем программу, описанную в предыдущем разделе.

Добавим в блок инструкций #include инструкцию, подключающую библиотеку pgraph. Теперь блок будет выглядеть так:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "pgraph.h"

При построении графика плотности распределения вероятностей понадобится константа-макрос, равная приближённому значению 12π. Внесём соответствующее дополнение в блок инструкций #define:

#define sqrt_2 1.414213562373       //sqrt(2)
#define sqrt_pi_2 1.253314137316    //sqrt(pi / 2)
#define koef_exp 0.398942280401     //1 / sqrt (2 * pi)
#define eps 1e-5                    //погрешность

Глобальную константу it инициализируем сразу же в момент объявления:

int it = 0;

Определение функция calc() нам больше не нужно. Заменим его определением функции get_norm_rand_val().

Единственное изменение функции f_1() заключается в удалении из её кода 17-й строки, в которой инициализируется нулём переменная it. Теперь it после формирования выборки будет содержать общее количество итераций, потребовавшихся для генерации сразу всех чисел.

Старую версию функции main() удаляем и на её место помещаем новую.

Перед тем, как привести её код, объясню принцип построения гистограммы.

Как мы выяснили, все значения, возвращаемые функцией get_norm_rand_val(), укладываются в отрезок [−4,3; 4,3]. Но нам будет удобнее работать с промежутком [−4,5; 4,5). Разобьём его на следующие 18 полуинтервалов одинаковой длины:

[−4,5; −4), [−4; −3,5), …, [4; 4,5).

Далее 10000 раз вызываем функцию get_norm_rand_val() и подсчитываем, сколько чисел, возвращаемых ею, попадает в каждый из 18-ти полуинтервалов. Для таких подсчётов создаём "массив счётчиков", состоящий из 18-ти элементов, каждый из которых соответствует "своему" полуинтервалу.

Изначально все элементы массива счётчиков инициализированы нулями. После получения каждого псевдослучайного числа выясняем, в какой полуинтервал он попадает, и увеличиваем на единицу значение соответствующего счётчика.

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

Далее наносим на холст для рисования координатные оси и для каждого полуинтервала вида [pq) строим прямоугольный столбик, основанием которого служит отрезок [pq]. Высота столбика равна отношению доли чисел, попавших в полуинтервал [pq) (т. е. частному от деления количества чисел, попавших в полуинтервал, на их общее количество, т. е. на 10000; такое частное называется относительной частотой) к длине полуинтервала, равной 0,5.

В результате получаем гистограмму, состоящую из 18 столбиков, сумма площадей которых равна 1. Единице равна и площадь под графиком плотности распределения вероятностей, который мы тоже построим, "наложив" его на гистограмму.

А вот и код новой версии функции main(), в которой реализованы все описанные выше действия:

int main()
{
    int arr[18];        //Массив счётчиков попаданий в промежутки
    for (int i = 0; i < 18; i++) //Инициализация элементов массива
        arr[i] = 0;              //счётчиков нулями
    //Генерация псевдослучайных чисел и подсчёт количеств попаданий в полуинтервалы
    for (int i = 0; i < 10000; i++)
        arr[(int) (9 + 2 * get_norm_rand_val())]++;
    //Создание холста для рисования
    image *img = create_image(581, 420);  //размер холста - 581x420 пикселей
    //Рисование координатных осей (начало координат - в точке (290; 9))
    line(img, 0, 9, 580, 9);        //Ось абсцисс (единица по оси абсцисс = 60 пикселей)
    line(img, 290, 0, 290, 419);    //Ось ординат (единица по оси ординат = 1000 пикселей)
    //Построение гистограммы
    set_cur_col(img, BLUE);            //Цвет - синий
    set_cur_pnt(img, 20, 9);           //Установка начальной точки
    for (int i = 0; i < 18; i++)
    {
        int h = round((0.2 * arr[i]));        //Высота столбика в пикселях
        line_to(img, img->cur_pnt.x, img->cur_pnt.y + h);    //Левая граница столбика
        line_to(img, img->cur_pnt.x + 30, img->cur_pnt.y);   //Верхняя граница столбика
        line_to(img, img->cur_pnt.x, img->cur_pnt.y - h);    //Правая граница столбика
    }
    //Построение графика плотности
    set_cur_col(img, RED);        //Цвет - красный
    set_cur_pnt(img, 20, 9);      //Установка начальной точки
    for (int i = 21; i < 561; i++)    //i - абсцисса текущей точки в сист. коорд. холста
    {
        double x = (i - 290) / 60.0;            //Абсцисса текущей точки
        double y = koef_exp * exp(-x * x / 2);  //Ордината текущей точки
        int j = round((y * 1000) + 9);    //Ордината текущей точки в сист. коорд. холста
        line_to(img, i, j);               //Проведение линии из предыдущей точки в текущую
    }
    //Сохранение изображения в файл
    save_to_file(img, "result.bmp");
    //Удаление холста
    free(img);
    //Печать общего количества итераций
    printf("Всего итераций: %d", it);
    return 0;
}

Думаю, дальнейшие пояснения излишни, поскольку комментарии приведены в самом коде.

Компиляция и запуск программы приводит к выводу на консоль следующего текста:

Всего итераций: 33262

Итак, среднее количество итераций, приходящееся на один вызов функции get_norm_rand_val(), — 3,3262, т. е. меньше четырёх. Это очень неплохо!

Ещё одним результатом выполнения программы является появление в корневой директории исполняемого файла графического файла result.bmp, содержащего следующее изображение:

Гистограмма и график плотности распределения стандартной нормальной случайной величины

Гистограмма и график плотности распределения вероятностей

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

Окончательный вид функций, генерирующих нормально распределённые псевдослучайные числа

Ниже приведён код, позволяющий генерировать нормально распределённые псевдослучайные числа, "очищенный" от переменных, функций и других объектов, предназначенных для тестирования.

#include <math.h>

#define sqrt_2 1.414213562373         //sqrt(2)
#define sqrt_pi_2 1.253314137316      //sqrt(pi / 2)
#define eps 1e-5

const double val[] = {0.0,            //f_1(0.5)
                      0.253347,       //f_1(0.6)
                      0.524401,       //f_1(0.7)
                      0.841621,       //f_1(0.8)
                      1.281552,       //f_1(0.9)
                      1.880794};      //f_1(0.97)

double f_1(double x)
{
    if (x == 0.5)
        return 0.0;
    if (x < eps)
        return -4.3;
    if (x > 1 - eps)
        return 4.3;
    int sign = 1;
    if (x < 0.5)
    {
        x = 1 - x;
        sign = -1;
    }
    double xk1 = val[(int) (round(10 *x) - 5)], xk;
    double b = 1 - 2 * x;
    do
    {
        xk = xk1;
        double t = xk / sqrt_2;
        xk1 = xk - sqrt_pi_2 * exp(t * t) * (erf(t) + b);
    }
    while (fabs(xk - xk1) >= eps);
    return sign == 1 ? xk1 : -xk1;
}

double get_norm_rand_val()
{
    return f_1(rand() / ((double) RAND_MAX));
}

Генерация псевдослучайных чисел, распределённых по нормальному закону с произвольными параметрами m и σ

Если X — случайная величина, распределённая по стандартному нормальному закону, то несложно показать, что случайная величина Y, выражающаяся через X следующим образом:

Y = m + σ X,

распределена по нормальному закону с параметрами m и σ. Таким образом, если, например, в некоторой программе значения параметров m и σ содержатся в переменных m и sigma соответственно, то для генерации псевдослучайных чисел, нормально распределённых с параметрами m и σ, можно многократно осуществлять вызов функции get_norm_rand_val() внутри следующего выражения:

... (m + sigma * get_norm_rand_val()) ...;

Здесь многоточиями обозначен некоторый код, каким-либо образом обрабатывающий значение выражения.

Заключение

Итак, поставленная задача успешно решена, причём её решение оказалось получить достаточно несложно. Теперь связку функций f_1() и get_norm_rand_val() можно использовать везде, где требуется генерировать псевдослучайные числа, распределённые по нормальному закону.

Скачать исходный код