Форсайт Дж., Малькольм М., Моулер К. - Машинные методы математических вычислений (1040536), страница 45
Текст из файла (страница 45)
вая при а-м обращении к ней - 1'„„=- а)'„+с(гпо<) т), гг и 1. Эти числа обращаются в числа с плавающей точкой из интервала 1О,1), которые и будут значениями Г)КА)хР. Полученное значение 1'„„возвращается посредством параметра 1 г' и прк последующем вызове должно использоваться как фактический параметр. При первом обращении к ОКА!НР параметру 1У нужно присвоить в качестве начального произвольное целочисленное значение. Значения т, а и с вычисляются автоматически при первом входе в программу. Главное допущение, сделанное здесь, заключаегся в том, что машина использует двоччиое представление целого числа, а умножение выполняется по модулю т, где т — степень двойки.
Это предположение упрощает вычисления. Г)КА<<Р находит значение пг)2, проверяя псслсдовательные сте- ') ТО ЕСТЬ «УИИ<ЕР<аЛЬОЫй Дат.ИХ СЛУт йЬЫХ ИСЕЛ» И <Дгтю<к РаеИОЫЕРИО рисорелелеиьых слух< и< ь х «се<» — Прил . «рсв. го. ВЫРАБОТКА случайных чисел 2б4 пени двойки до тех пор, пока умножение на 2 уже больше пе дает приращения величины, Предполагается также, что сложение целых чисел выполняется по модулю т либо сохраняются по крайней мере !он,(т) значимых разрядов. Значения а и с вычисляготся в соответствии с указанными выше рекомендациями Кнута. В подпрограмме а соответствует параметр 1А, а с — параметр 1С.
Случайный шаблон разрядов для а достигается обращением РАТА!ч' (!.РО), дающим значение с двойной точностью числа и/4, которое на двоичной машине совпадает со сдвинутым шаблоном для и. Деление на 8.РО и умножение на ог/2 также, надо полагать, не изменят этого шаблона. Наконец, значение с двойной точностью преобразуется в целое число, умножается на 8 и увеличивается на 5, чтобы обеспечить условие а(шог) 8) — 5. Полученное значение а равно примерно (нт/8)пжгигг2. Зтим удовлетворяются ограничения в виде неравенств.
Значение с вычисляется прямо пз определения 3. Некоторые трансляторы с ФОРТРАНа не обращают константы типа 8.РО в точное представление с плавающей точкой, но это едва ли будет иметь серьезные последствия. Как доказано в теореме А (Кнут (1969), стр. 29), последовательность (1'„) обязательно имеет период максимальной длины рт. Однако наименее значимые двоичные разряды )т„будут «не очень случайнымн». Когда К„преобразуются в числа с плаваюшей точкой, то эти наименее значимые разряды обычно не важны.
Чтобы вычислять случайные целые числа между 0 и К вЂ” 1, нужно использовать фу нкцикт ! Р)Х(К»ЩА)т)Р(1У)). Для машин серии 360 фирма 1ВМ вместе со своей фортранной системой поставляет 85Р') — пакет подпрограмм для научных расчетов. Эти подпрограммы могут автоматически вызываться из трансляторов с ФОРТРАНа, В пакет входит и датчик случайных чисел под названием КА)т)Р().
Л4ье наспгоягггвльно призываем вос исггольэованть РКАХР, но не асггольэовать КАУР() по причинам, которые сейчас будут объяснены '). Прн гснериронанин )гп КА)т!РП полагает а= 65539 н с=О. Период полученной последовательности вполне удовлетворителсн — 2»", т. е. четверть всей области целых чисел. (Период последовательности, вырабатываемой подпрограммой БАКАИР на 1ВМ 360, ранен максимально возможному значению 2".) Трудности связаны с тем обстоятельством, что 65 539=2"+3. Возможно, это число было выбрано для ускорения умножения, но оно приводит к катастрофическим свойствам последовательности.
»1 То есгь «сгепнйс ноьгоннпе распахе,— Прим, перев. ») Эпп ~ ризыв относится, равумеегся, к американским читателям книги— Прим, перев. !В!, РАВНОМЕРНО РАСПРЕДЕЛЕННЫЕ ЧИСЛА 2бз Следующие выкладки проводятся по модулю 2"; =- (2'"'+ 6. 2" + 9) х„=- 16 (2" + 3) — 9! хА =- бх,„— 9хА. Итак хх,=-бхь„— 9хА для всех л. Как результат этого в последовательности, вырабатываемой подпрограммой (ААХР(), наблюдается крайне высокая корреляция между тремя подряд идущими случайными целыми числами. Например, каждый раз, когда х„близко к нулю, хь„близко (по модулю 2!') к бх>, Мы проверили от 10 000 до !00 ООО троек чисел (х>„ха+„ х„„), генерируемых обеими подпрограммами: НЛХ!Р13 и ()КАХА!Р. Эти тройки можно раскладывать по 1000 ящикам, соответственно тому, в какую из десяти равных частей интервала 10,2" — 1) попадет каждая из переменных.
(Это равносильно тому, что для проверки А)НЛХ!Р интервал (О,!) разделили на десять равных частей.) Еслг! числа, вырабатываемые датчиком, действительно случайные, то каждая тройка должна иметь одинаковые шансы попасть в любой из !000 ящиков. К распределениютроек поящикам можно применить критерий у! (Кнут (!969)). Результаты показали, что подпрограмма РКАХР вполне удовлетворительна, в то время как (хАХ! РР совершенно дискредитировала себя в качестве датчика случайных чисел. Мы полагаем, что КАХР1) настолько плоха, что во многих случаях моделирование, основанное на ней, вероятно, имело определенный режим, особенно если использовались тройки чисел.
Поскольку ЙАХИЗ почти неизменно применяется пользователями Системы 360, то мы серьезно сомневаемся в результатах большого числа моделирований. Марсалья (1968) доказал, что все датчика случайных чисел, использующие рекуррентные соотношения, в определенной степени страдают корреляцией между последовательными числами, однако для хорошо составленных датчиков, таких, как ()1(ЛХР, этот эффект гораздо меньше, чем для многих других. Было изобретено много других тестов для датчиков случайных чисел, но обычно самая лучшая проверка датчика состоит в применении его к модельной задаче (из вашей области ппнложений), ответ которой вам известен.
В настоящее время растет интерес к методам построения последовательностей чисел х„из 10,1), которые не выглядят случайными, однако оказываются хорошо приспособленными для конкретного применения. Например, существуют последовательности хь, которые равномерно распределены на интервале 10,11) с гораздо большей регулярностью, чем подлинно случайные чис- ~о вьшавогка случайных чисел 266 ла. Онн могут быть полезны для интегрирования по методу Монте-Карло вследствие меньшей дисперсии оценок.
Если скорость выработки равномерно распределенных случайных чисел жизненно важна для вашей программы, вам нужно изучить специализированные программы для вашего класса машин, публикуемые в текущей литературе. Например, имеется программа на ФОРТРАНе (Марсалья, Брэй (1968)) и программа для машин серии 360, написанная в машинном коде (Серафин (1969)). Время выполнения последней на машине 360767 составляет 11.9 микросекунд плюс 4.2 микросекунды на вызов подпрограммы. Для подпрограммы 13КА[т[О на ! ВМ 360/67 требуется около 90 микросекунд (!аОКТРА[т[ (Н), (ОРТ=-2)).
40.2. Подпрограмма Ц[аА['[Р Следующая программа иллюстрирует тривиальное применение [)КА[4Р. Печатается десять строк; каждая из них содержит в случайном порядке слова )!еабз и (а!Ы. С Иллюстрирующая программа для ОКЛМР С 1У=О РО !О 1 —. 1, 1О 1Р(()КАМР(1Ъ') Л.Т. О. 6) ОО ТО 5 ЖК!теу(6, Й ОО ТО 1О 6 !УК П Е(6, 2) 10 СОМТ1ЬЦЕ ГОКМЛТ(ОН НЕАР5) 2 РОКМАТ(ОН ТАР.5) 5ТОР ЕМР Другое начальное значение для переменной [Ъ' приведет к другой последовательности случайных чисел. Более того, даже при одном и том же начальном значении для 1У иа машинах с различным размером слова будут получены разные последовательности.
КЕЛЬ ГЦМСТ)ОХ РКЛКР(1У) 1мтебеК 1У С С ЦКЛКР— ЭТО ДАТЧИК РАВНОМЕРНО РАСПРЕДЕЛЕННЫ С СЛУЧАЙНЫХ ЧИСЕЛ, ОСНОВАННЫЙ НА ТЕОРИИ И ПРЕДЛ С ЖЕНИЯХ, СОДЕРЖАШИХСЯ В КНИГЕ КНУТ (!ЧОЗ). ТОМ С ПЕРЕД ПЕРВЫМ ОБРАЩЕНИЕМ К ЕКАМР ЦЕЛОЙ ПЕРЕМЕ Г нОИ !у слеДУет присВОить пРОизВОльнОе ЦелОчисле С НОЕ НАЧАЛЬНОЕ ЗНАЧЕНИЕ. ВЫЗЫВАЮК(АЯ ПРОГРЛММ С НЕ ДОЛЖНЛ ИЗМЕНЯТЬ ЗНАЧЕНИЕ (У МЕЖДУ ПОСЛЕДОВ Г ТЕЛЬНЫМИ ВЫЗОВЛМИ. ЗНЛс!ЕНИЯ ФУНКЦИИ ЦКАМР ЯВЛ С 1ОТСЯ ЧИСЛЛМИ ИЗ ИНТЕРВАЛА (О, 1). С Ш БЫРЛБОТКЛ СЛУЧЛЙНЫХ ЧИСЕЛ 10.3. Выбор иа других распределений Часто выработка случайного числа из [0,1) является просто средством принятия некоторого случайного решения. В качестве очень простого примера предположим, что в некой программе моделирования желательно при случайном выборе попадать на какую-то ветвь одну десятую часть всего времени.
Этого легко добиться, выбирая случайное число у на [0,1) и беря данную ветвь, если у<0.1. Более быстрый способ — выбирать случайное целое число в [О, 2"' — 1) и брать эту ветвь при х~2"/1О. Таким же образом из равномерно распределенных случайных целых чисел можно получить любое конечное распределение с различными весами.
Гораздо более сложно вычисление случайного числа из неравномерного непрерывного распределения. Если для функции распределения Е(х) не имеется каких-либо специальных приемов, то общий метод заключается в следующем: берут случайное число, равномерно распределенное на [О,1), а затем, пользуясь обратной функцией, образуют у=-Г-'(х). Трудной проблемой при этом является достаточно быстрое и достаточно точное формирование Е-'. Некоторые часто используемые и важные распределения уже довольно полно исследованы. Вероятно, наиболее известным является нормальное распределение й[ (0,1) с нулевым средним значением и дисперсией 1. Оперировать с обратным отображением Г ' сравнительно трудно и медленно.
Распространенный подход таков: берут 12 равномерно распределенных на [0,1) чисел уь затем по формуле х;=2у; — 1 получают 12 случайных чисел, равномерно распределенных на [ — 1,1), и складывают их. Результат будет довольно хорошим приближением к выбору из распределения Лl (0,1).
Метод является относительно медленным. Намного более изобретательный метод принадлежит Боксу и Мюллеру. Приводим его в более поздней модификации Джеймса Белла из Стэнфорда. 1. Образовать два случайных числа (1, и ([,, равномерно распределенных на [0,1) (например, пользуясь подпрограмьюй О[х АХВ). 2. Образовать [У,==2(1,— 1 и Р.,— 2(/,— 1; эти числа равно-1 мерно распределены на [ — 1,1).
3. Образовать 5 — у[+У,,'. Если 5'- 1, отбросить [У, и Г, и 1 вернуться к шагу 1. (Здесь мы теряем 22% эффективности.)~ Если 5(1, то (Ро Ре) — случайная точка, равномерно распре- деленная в единичном круге. 4. Образовать эГ ~!п~ 1( [у эà — Э!Я5 УПРАЖНЕНИЯ 269 где !и — натуральный логарифм (на зтот шаг затрачивается ббльшая часть времени).