Моделирование нормально распределенных случайных величин

Лекция 25. Моделирование нормально распределенных случайных величин Нормальный закон распределения встречается в природе весьма часто, поэтому для него разработаны отдельные эффективные методы моделирования. Формула распределения вероятности значений случайной величины x по нормальному закону имеет вид: Как видно, нормальное распределение имеет два параметра: математическое ожидание mx и среднеквадратичное отклонение σx величины x от этого математического ожидания. x — случайная величина; y(x) — вероятность принятия случайной величиной значения x; mx — математическое ожидание; σx — среднее квадратичное отклонение. Нормализованным нормальным распределением называется такое нормальное распределение, у которого mx = 0 и σx = 1. Из нормализованного распределения можно получить любое другое нормальное распределение с заданными mx и σx по формуле: z = mx + x · σx. Рассматривая последнюю формулу, вспомните формулы компьютерной графики: операция масштабирования выражается в математической модели через умножение (это соответствует изменению разброса величины, растягиванию геометрического образа), операция смещения выражается через сложение (это соответствует изменению значения наиболее вероятной величины, смещению геометрического образа). Функция нормального распределения имеет вид колокола. На рис. 25.1 показано нормализованное нормальное распределение. Рис. 25.1. Графический вид нормального закона распределения случайной величины х с параметрами mx = 0 и σx = 1 (распределение нормализовано) 1 График на рис. 25.1 показывает, что в области –σ < x < σ на графике сосредоточено 68% площади распределения, в области –2σ < x < 2σ на графике сосредоточено 95.4% площади распределения, в области –3σ < x < 3σ на графике сосредоточено 99.7% площади распределения («правило трех сигм»). Вспомните, пожалуйста, рис. 2.7 из лекции 02. Пример. По нормальному распределению распределен рост людей, находящихся одновременно в большой аудитории. А именно: достаточно мало людей очень большого роста, и столь же мала вероятность встретить людей очень малого роста. В основном, легче встретить людей среднего роста — и вероятность этого велика. Например, средний рост людей составляет, в основном, 170 см, то есть mx = 170. Известно также, что σx = 20. На рис. 25.1 показано, что доля людей с ростом от 150 до 190 (170 – 20 < 170 < 170 + 20) составляет в обществе 68%. Доля людей от 130 см до 210 см (170 – 2 · 20 < 170 < 170 + 2 · 20) составляет в обществе 95.4%. Доля людей от 110 см до 230 (170 – 3 · 20 < 170 < 170 + 3 · 20) составляет в обществе 99.7%. Например, вероятность того, что человек окажется ростом меньше 110 см или больше 230 см составляет всего 3 человека на 1000. Свойства нормального распределения Изменение параметра нормального распределения mx приводит к сдвигу кривой по оси x (см.рис. 25.2). Рис. 25.2. Влияние параметра «математическое ожидание» на вид закона нормального распределения случайной величины х Изменение параметра нормального распределения σx приводит к масштабированию формы по оси x (напоминаем, в любом случае всегда площадь под кривой плотности вероятности неизменна и равна 1). Чем более не случаен процесс, тем меньше его среднеквадратичное отклонение, тем уже и выше колокол на графике. Изменение параметра нормального распределения σx приводит к масштабированию формы (см. рис. 25.3) по оси x (напоминаем, в любом случае всегда площадь под кривой плотности вероятности неизменна и равна 1). 2 Рис. 25.3. Влияние параметра «среднеквадратичное отклонение» на вид закона нормального распределения случайной величины х Чем более не случаен процесс, тем меньше его среднеквадратичное отклонение, тем уже и выше колокол на графике. Действительно, разброс случайности относительно математического ожидания становится все более минимальным. В пределе детерминированный процесс имеет вид, показанный на рис. 25.4. Рис. 25.4. Вид закона нормального распределения вероятности при переходе его к детерминированному случаю в пределе (σx = 0). Случайное событие становится детерминированным: x = mx ± 0 (разброса нет) Изучать детерминированные процессы проще. Чем больше величина σx, тем менее закономерно поведение изучаемого объекта, так как возможны любые значения характеризующих его параметров, разброс величин относительно средней ожидаемой увеличивается. Прогнозирование и управление поведением объекта в этом случае затрудняется. 3 Рассмотрим вид интегральной кривой плотности распределения случайной величины, распределенной по нормального закону. Вид ее приведен на рис. 25.5. F — интегральная функция Лапласа. Смысл интегральной функции — вероятность того, что случайная величина примет значения из диапазона от –∞ до x. Например, запись F(170) = 0.5 для нашего примера означает: вероятность того, что случайно выбранный из аудитории человек будет ростом не выше 170 см, составляет 0.5 (то есть каждый второй). Рис. 25.5. Вид интегральной функции Лапласа F(x) Данная функция задана интегралом от плотности вероятности нормального распределения: К сожалению, этот интеграл не берется в общем виде, поэтому функция Лапласа задана в видетаблицы для mx = 0 и σx = 1. Поскольку функция Лапласа симметрична относительно точки(x = 0, y = 0.5) (как и функция самого нормального распределения), F(–x) = 1 – F(x), то в таблице содержится только одна из ее симметричных частей. Если задается интервал интегрирования функции Лапласа [a; b], то: Вероятность попадания X в интервал, симметричный относительно mx: Например, для правила «трех сигм»: P(|x – mx| < 3σ) = 2 · F(3) – 1 = 2 · 0.9987 – 1 = 0.9973 (как ранее мы и указывали). Число F(3) = 0.9987 взято из таблицы Лапласа. Пример. Найти вероятность изготовления детали с ошибкой в ее размерах не более 15 мм, если известно, что изготовление детали с ошибкой распределено по нормальному закону m = 0 и σ = 10 мм. 4 P(|x| < 15) = P(–15 < x < 15) = F((15 – 0)/10) – F((–15 – 0)/10) = F(1.5) – F(–1.5) = F(1.5) – (1 – – F(1.5)) = 2 · F(1.5) – 1 = 2 · 0.9332 – 1 = 0.8664. То есть 8664 деталей из 10000 будут иметь ошибку в размерах не более 15 мм. Важным свойством закона является то, что нормальное распределение является пределом для различного вида распределений, вытекающее из центральной предельной теоремы: «Для большого числа N случайных величин X с любым законом распределения их сумма есть случайное число с нормальным законом распределения». где a — математическое ожидание в законе распределения случайной величины X; σ — среднеквадратичное отклонение в законе распределения случайной величины X; N — количество случайных чисел. Вспомним опыт «Доска Гальтона» из физики (см. рис. 25.6). Рис. 25.6. Доска Гальтона. Шарики, падающие сверху в сосуд случайно распределяются в нем в согласии с нормальным законом распределения Доска разделена на секции; в верхней части доски находятся особым образом расположенные стержни, ударяясь о которые, множество падающих сверху под влиянием силы тяжести шариков, испытывая соударения также и между собой, меняют свою траекторию полета. В результате в различные секции попадает различное количество шариков. Если дождаться окончания опыта, то в каждой секции будет определенное количество шариков, конечно, каждый раз разное, так как процесс их столкновения случайный. Но интересно то, что распределение шариков по секциям будет образовывать нормальный закон распределения. Вроде бы доска не меняется, шарики падают одни и те же, и тем не менее, во-первых, форма распределения слегка колеблется 5 (случайность), и разные шарики попадают в разные секции, во-вторых, на макроуровне, где проявляется организация шариков как совокупности, всегда получается нормальный закон распределения (закономерность). Для пытливого ума опыт порождает очень интересные вопросы. Почему малые отклонения в поведении элементов системы (шариков) ведут к большим разбросам в их координатах? Почему случайности не компенсируются? Насколько они не компенсируются? Обозначим через D расстояние между стержнями доски Гальтона, r — диаметр стержня (см.рис. 25.7, а, б). Очевидно, что отклонение β шарика от отвесной траектории при малом угле α(изменение направления вследствие удара о стержень) можно посчитать как tg(α1) = β1/r или α1 ≈ β1/r(при малых углах α1). То есть, далее шарик летит под углом α1 к траектории свободного падения на расстояние D, пока не встретит новый стержень. За время перелета по прямой от стержня к стержню шарик отклонится от отвесной вертикальной линии на величину β2 = α1 · D. Далее снова, испытав соударение со стержнями, шарик отклонится на угол α2 ≈ β2/r и β3 ≈ α2 · D = β1 · (D/r)2 и так далее. Рис. 25.7. Взаимодействие шариков со стержнями доски (детализация) Если шарик испытает n соударений, то он отклонится от отвесной траектории на расстояниеβn ≈ β1 · (D/r)n – 1 (конечно, если ему «повезет» и он будет отклоняться все время в одну сторону). Теперь посмотрим, сколько это означает в цифрах, например, D = 20r, пусть стержней, от которых отклонился шарик в полете (соударений), n = 10, β1 = 1 нм (расстояние меньше атома!!), то тогда легко вычислить, что β10 = 10.24 км. То есть, видно, что сверхмалые отклонения (1 нм и число соударенийn = 10), фактически случайности, приводят к макроэффектам. На самом деле такой разброс наступает гораздо быстрее, так как шарики сталкиваются еще и между собой. Кстати, удаление какого-то шарика на 10 км — это тоже невероятное событие, так как вероятность, что шарику «повезет» 20 раз подряд, составляет крайне маловероятное число P = 0.520 = 10– 6. Рассмотрим еще один эксперимент «Распределение Максвелла». В сосуде находится 40 синих и 40 красных шаров, каждые в своей половинке сосуда. Половинки сосуда разделены перегородкой с отверстием, через которое могут проникать шары из одной части сосуда в другую (см. рис. 25.8). Количество шаров, как синих, так и красных, в каждой из частей сосуда подсчитываются. 6 Рис. 25.8. Схема эксперимента, демонстрирующего диффузию (появление хаоса) в сложных системах Шары имитируют броуновское движение молекул, соударяясь друг с другом упруго, обмениваясь при ударе энергиями, импульсами, направлениями скорости и также взаимодействуя со стенками сосуда. При смоделированном в опыте упругом ударе усредняется средняя кинетическая энергия встречающихся частиц. Несмотря на то, и в этом суть диффузии, что шары находились каждый цвет строго в своей половине, со временем примерно половина (50%) красных шаров окажется в первой половине сосуда, а половина (50%) — во второй его части, то же касается и синих шаров (50% : 50%). Употребленное здесь слово «примерно» означает, что время от времени совокупность шаров будет делиться не ровно на половины, а в пропорции (50% ± h : 50% ± h), так как случайное количество шаров будет переходить в произвольные моменты времени из одной части сосуда в другую, и наоборот. Обнаруженная в наших опытах относительная флуктуация средних энергий достаточно велика, порядка 20%. Несмотря на то, что шары каждого цвета находились строго в своей половине (в этом и состоит суть диффузии), со временем примерно половина (50%) красных шаров окажется в первой половине сосуда, а половина (50%) — во второй его части, то же касается и синих шаров (50%:50%) (см. рис. 25.9). Рис. 25.9. Типичный вид поведения характеристик в системах с диффузией. На рисунке показано изменение количества синих и красных молекул в одной из половинок сосуда со временем Особенность 1. Можете продолжать наблюдение сколь угодно долго, но при этом никогда первоначального распределения этих шаров достигнуто не будет. И, наоборот, при любом 7 начальном порядке шаров система будет переходить к одному и тому же хаотическому состоянию. Почему? Объяснение этого надо искать в том, что время необратимо. Операция усреднения также необратима. Если я вам сообщу, что сумма двух чисел равна 8, то вряд ли вы сможете сказать, какие я выбрал для этого слагаемые. И множество вариантов слагаемых ведут к одной и той же сумме. Хаотическая система не помнит свое прошлое. А именно: состояние такой системы по прошествии времени не зависит от того, как и из какого состояния система пришла в хаотическое. Из любого начального состояния наша система довольно быстро переходит в это хаотическое состояние, и находится в нем далее сколько угодно долго. То есть хаотическое состояние — устойчивое состояние. Интересно, что в этом устойчивом состоянии можно отыскать такой ряд параметров, который будет неизменен со временем, такие параметры описывают систему в целом и называютсямакропараметрами. Например, в нашем примере устойчивым макропараметром является распределение количества молекул n по модулю скорости v, которое демонстрирует закон Максвеллаn(v) = r · v2 · e–c · v2/T (r, c — коэффициенты, T — температура совокупности молекул). Особенность 2. Если разделить сосуд на большее число частей, то шары распределятся между ними примерно в равных долях. Особенность 3. Слово «примерно» пополам, говорит нам о том, что в ответе 50% : 50% всегда есть разброс, дисперсия, флуктуация. Этот разброс уменьшается, если увеличить число молекул в сосуде. Особенность 4. Если запустить модель снова и попробовать шары сначала упорядочить, например, установить их ровно рядами в одной половине сосуда, или вверху, или в одном из углов, то очень быстро шары «примерно» равномерно распределятся по сосуду, и наступит хаос, то есть порядок исчезнет и, что интересно, никогда снова не наступит. Откуда возникает хаос, как система переходит от порядка к хаосу? Причина хаоса, как мы уже ранее показали в опыте «Доска Гальтона», появляется в момент небольших отклонений шаров (пусть, например, 10–9 м) от идеальной теоретической траектории, но это небольшое отклонение на пролете расстояния D приводит к существенному отклонению траектории от идеала и столкновению шаров, что резко меняет траектории соседей. Мгновенно появляется лавина соударений, и хаос наступает очень быстро. Особенность 5. Ранее, считалось, что системы с малым числом элементов не способны проявлять случайность, хаотическое поведение. Преобладало рассуждение физиков о газе из миллиардов молекул. Проведите этот же эксперимент с пятью шарами, потом с тремя и, наконец, одним шаром. Распределение Максвелла в результате эксперимента проявляется достаточно упрямо, правда, для этого требуется ждать в несколько раз дольше. (Исключение составит опыт с одним шаром!) Опыт показывает, что усреднение по ансамблю молекул и по времени проводит к одинаковому результату. Так что хаос может возникнуть и на небольших множествах. Итак, и опыт «Доска Гальтона» и экспериментальное снятие распределения Максвелла приводят нас эмпирически к тому, что нормальное распределение широко распространено в окружающем нас мире. Суммирование случайных величин, сумма действий, совместное наложение эффектов на микроуровне часто приводит, усреднившись, к появлению нормального распределения на макроуровне. Рассмотрим вопрос имитации случайных величин, заданных нормальным законом распределения. Табличный метод генерации нормально распределенных чисел 8 Для этого нормальное число можно взять из справочника в таблице функции Лапласа и получить случайное число по методу взятия обратной функции (см. лекцию 24): x = F–1(r), где F — интегральная функция Лапласа. Технически это означает, что надо разыграть случайное равномерно распределенное число r из интервала [0; 1] стандартным ГСЧ (см. таблицу абсолютно случайных проверенных чисел), найти равное ему число в таблице значений функции Лапласа в столбце F и по строке определить случайную величину x, соответствующую этому числу. Недостатком метода является необходимость хранения в памяти компьютера всей таблицы чисел функции Лапласа. Метод генерации нормально распределенных чисел, использующий центральную предельную теорему Общая идея метода следующая: требуется сложить случайные числа с любым законом распределения, нормализовать их и перевести в нужный диапазон нормального распределения. Допустим, что нам надо в целях имитации получить ряд случайных чисел x, распределенных по нормальному закону с заданными математическим ожиданием mx и среднеквадратичным отклонениемσx. 1. Сложим n случайных чисел, используя стандартный ГСЧ: Согласно ЦПТ числа V образуют ряд значений, распределенный по нормальному закону. Эти числа тем лучше описывают нормальный закон, чем больше параметр n. На практике n берут равными 6 или 12. Заметим, что закон распределения чисел V имеет математическое ожиданиеmV = n/2, σV = sqrt(n/12). Поэтому он является смещенным относительно заданного произвольного. 2. С помощью формулы z = (V – mV)/σV нормализуем этот ряд. Получим нормализованный закон нормального распределения чисел Z. То есть mz = 0, σz = 1. 3. Формулой (сдвиг на mx и масштабирование на σx) преобразуем ряд Z в ряд x: x = z · σx + mx. Пример. Смоделировать поток заготовок для обработки их на станке. Известно, что длина заготовки колеблется случайным образом. Средняя длина заготовки составляет 35 см, а среднеквадратичное отклонение реальной длины от средней составляет 10 см. То есть по условиям задачи mx = 35, σx = 10. Тогда значение случайной величины будет рассчитываться по формуле:V = r1 + r2 + r3 + r4 + r5 + r6, где r — случайные числа из ГСЧрр [0; 1], n = 6. X = σx · (sqrt(12/n) · (V – n/2)) + mx = 10 · sqrt(2) · (V – 3) + 35 или X = 10 · sqrt(2) · ((r1 + r2 + r3 + r4 + r5 + r6) – 3) + 35. Метод Мюллера Совсем простым методом получения нормальных чисел является метод Мюллера, использующий формулы: Z = √(–2 · Ln(r1)) · cos(2π · r2), где r1 и r2 — случайные числа из ГСЧрр [0; 1]. Можно также воспользоваться аналогичной формулой Z = √(–2 · Ln(r1)) · sin(2π · r2), где r1 и r2 — случайные числа из ГСЧрр [0; 1]. 9 Пример. Материал поступает в цех один раз в сутки по 10 штук сразу. Расход материала из цеха случайный по нормальному закону с математическим ожиданием m = 10 и среднеквадратичным отклонением σ = 3.5. Вычислить вероятность дефицита на складе при запасе материала в начальный момент времени 20 штук. При реализации в среде моделирования Stratum решение задачи будет выглядеть следующим образом. z := sqrt(–2 · ln(rnd)) · cos(2 · PI · rnd) x := σ · z + m запас := запас + 10 запас := (запас – x) · ed(запас – x) дефицит := дефицит + not(ed(запас – x)) N := N + 1 P := дефицит/N stop(N > k) z — нормальное нормализованное случайное число; x — нормальное число, ежедневный расход материалов; запас — состояние склада: начало дня, моделирование прихода; запас — состояние склада: конец дня, моделирование расхода; дефицит — счетчик дней, в течение которых наблюдался дефицит; N — количество дней; P — вероятность дефицита; k — моделирование в течение k дней. 10