Лабораторная работа № 1

Разработка генератора псевдослучайных чисел.

Определение длины периода и отрезка апериодичности.

проверка гипотезы о равномерности распределения случайной последовательности по критерию χ²

Общие положения

В задачах, решаемых методом Монте- Карло, используются случайные числа. Однако для решения задачи нам не важно, как получены эти числа, лишь бы они удовлетворяли тестам: были независимыми и одинаково распределенными. Одним из способов получения таких чисел является метод псевдослучайных чисел. Числа γ₁, γ₂,… γ_n , которые вычисляются по какой-либо заданной формуле и могут быть использованы вместо случайных чисел при решении некоторых задач, называются псевдослучайными числами. Достоинством этого метода является то, что существуют простые формулы для расчета этих чисел, причем на получение каждого числа затрачивается 3-5 команд ЭВМ, а программа расчета занимает в накопителе лишь несколько ячеек.

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

γ_n+1=Ф(γ_n), (1)

где начальное число γ₀ задано.

Например, функция y=Ф(x), изображенная на рис.1, не может породить хорошую последовательность псевдослучайных чисел γ₁, γ₂,… γ_n. Действительно, все пары чисел (x,y) лежат на кривой y=Ф(x), а равномерно распределенные пары должны равномерно заполнять квадрат 0<x<1, 0<y<1.

Следовательно, “хорошую” последовательность может породить только такая функция y=Ф(x), график которой весьма плотно заполняет единичный квадрат. Примером такой функции может служить функция y=Д(Ax) при больших значениях А, где Д(u) –дробная часть числа u. График функции γ_n+1=Д(12γ_n) изображен на рис.2.

Рассмотрим некоторые процедуры получения последовательностей псевдослучайных равномерно распределенных чисел.

y y y=Д(12x)

1 1

Y=Ф(x)

1 x 0 1 x

Рис.1 Рис.2

Первый алгоритм для получения псевдослучайных чисел был предложен Дж. Нейманом. Это метод середины квадратов (middle- square method). В этом методе число γ_n предполагается 2k- значным:

γ_n=0,a₁,a₂,…,a_2k.

Чтобы получить число γ_n+1, надо γ_n возвести в квадрат

γ_n ² =0,b₁,b₂,…,b_4k

и затем отобрать средние 2k цифр этого квадрата:

γ_n+1= 0,b_k+1,b_k+2,…,b_3k.

Нетрудно проверить, что этому методу соответствует функция

Ф(x)=Д(10^-2k Ц(10^3kx² ) ), (2)

или, что то же

Ф(x)=10^-2лЦ(10^2kД(10^kx²)),

где Ц(u)- целая часть числа u : u=Д(u)+Ц(u).

Однако в последовательностях, построенных этим методом, получается больше, чем нужно, малых чисел. А при некоторых γ₀ наблюдается вырождение последовательности.

Наибольшее распространение получил алгоритм, предложенный Д. Лемером, который называется методом вычетов (residual method), или методом сравнений (congruential method).

В этом методе Ф(х)=Д(gx), т.е.,

γ_n+1=Д(gγ_n), (3)

где g– большое целое число.

Если задать γ₀ в форме несократимой дроби γ₀ =m₀/M, где m₀ и M – целые числа и М взаимно просто с g, то все γ_n будут несократимыми дробями вида γ_n = m_n/M, где числители m_n определяются формулой

m_n+1=gm_n(mod M) (4)

Запись (4) означает, что m_n+1 равно остатку, полученному при делении gm_n на M ( или, другими словами, m_n+1- это наименьший положительный вычет gm_n по модулю M). При некоторых g, M, m₀ получаются удовлетворительные последовательности, при других- плохие.

Удовлетворительная последовательность псевдослучайных чисел получается , например, при g=5¹⁷, M=2⁴², m₀=1.

Еще один пример алгоритма типа (2)

γ_n+1=10^-kЦ(10^kД((1- γ_n)³10^k)), например, k=5.

Возможна также модификация метода вычетов путем прибавки к числителю некоторой константы

m_n+1=(gm_n+C) (mod M) .

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

Единственный недостаток метода – ограниченность количества псевдослучайных чисел, ибо, если последовательность чисел γ_0, γ_1,… γ_k… вычисляется на ЭВМ по формуле (1), то она обязательно периодическая. В самом деле, так как в коде любой ЭВМ можно записать лишь конечное число N чисел, заключенных между нулем и единицей, то рано или поздно какое-нибудь значение γ_L совпадет с одним из предыдущих значений γ_l . Тогда в силу (1)

γ_L+i= γ_l+i при i=1,2,… (5)

Пусть L- наименьшее число, удовлетворяющее (5) при некотором l (l<L). Множество чисел γ₀, γ₁,… γ_L-1 называется отрезком апериодичности последовательности (1), число L- длиной отрезка апериодичности, а P=L-l – длиной периода (рис.3)

Очевидно, что отрезок апериодичности состоит из различных чисел. И обычно для расчета не рекомендуется использовать больше, чем L чисел последовательности (1), т.к. это ведет к повторению испытаний в тех же условиях, что и раньше. Здесь увеличение числа реализаций не дает новых статистических результатов.

P P

0 1 2 … l l+1 … L-1 L… … i

L

Рис.3

Критерии проверки гипотез о законе распределения.

Случайные числа, полученные алгоритмическим способом, необходимо проверить на соответствие их предполагаемому закону распределения. Существует много критериев для такой проверки (критериев согласия):: критерий Колмогорова, критерий ω² и т.д. Каждый из них имеет свою область применения.

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

Рассмотрим критерий χ² (критерий Пирсона).

Критерий χ² (хи-квадрат) имеет очень широкое распространение в практике статистических исследований. Он определяется формулой

,(6)

где r- количество интервалов гистограммы, построенной по статистическому материалу;

m_i – количество реализаций, попавших в i–ый интервал;

n – общее количество реализаций, содержащихся в экспериментальных данных;

p_i – вероятность попадания случайной величины в i–ый интервал:

p_i=F(x_i+1) – F(x_i), (7)

где F(x_i+1) , F(x_i) – значение теоретической функции распределения соответственно на правой и левой границах i-го интервала гистограммы.

При расчетах принимают

F(x₁)=0,

F(x_r+1)=1.

Пирсон показал, что при увеличении объема выборки выборочное распределение величины χ² стремится к распределению, задаваемому плотностью

.

Это распределение называют распределением χ² с (r-1) степенями свободы.

Рекомендуется следующий порядок проверки гипотезы.

1. Результаты наблюдений случайной величины располагают в вариационный ряд:

x₁, x₂, …, x_n, ( х₁<x₂<…< x_n).

2. Размах вариационного ряда делят на r интервалов.

Величину r выбирают, пользуясь формулой Стаджерса

r= 1+3.3lg n,

или следующими рекомендациями:

при n≥ 1000 r=11…15,

при n≥ 400 r=10,

при n≥ 200 r=9.

При 100< n < 200 критерий применяют в исключительных случаях с числом интервалов r= 7…9.

Если интервалы выбраны одинаковыми, то ширина их равна .

3. Подсчитывают количество реализаций, попавших в каждый из интервалов.

4. Строят гистограмму.

5. По виду гистограммы формулируется гипотеза о виде закона распределения.

6. По статистическим данным вычисляются оценки параметров закона распределения.

7. Найденные в п.6 оценки принимаются в качестве параметров закона распределения случайной величины Х и вычисляются теоретические вероятности попадания случайной величины в каждый из интервалов гистограммы по формуле (7).

8. Определяется величина χ² по формуле (6).

9. По таблице распределения для вычисленного значения χ² и числа степеней свободы s=r-k-1, где k- количество параметров теоретического закона распределения (для экспоненциального к=1, для нормальноо и Вейбулла к=2), находится значение вероятности P того, что величина χ² примет значение, большее χ². Если эта вероятность большая, то считается, что нулевая гипотеза хорошо согласуется с результатами эксперимента, если маленькая, то не согласуется.

Таблица значений χ² в зависимости от P и s приведена ниже.