- •Случайные события
- •Стохастический мир
- •Случайные величины
- •Совместная и условная вероятности
- •Зависимость и независимость
- •Характеристическая функция
- •Модель аддитивного блуждания
- •Случайные процессы
- •Стохастические уравнения
- •Уравнение Ито
- •Остановка перед восхождением
- •Лемма Ито
- •Точные решения
- •Простые стохастические модели
- •Представление решений
- •Автокорреляция и спектр
- •Порождающий процесс Винера
- •Средние значения
- •Динамическое уравнение для средних
- •Процесс Феллера
- •Логистическое уравнение
- •Вероятности
- •Марковские плотности вероятности
- •Граничные условия
- •Стохастические интегралы
- •Площадь под траекторией Винера
- •Интегралы Ито
- •Интегрирование стохастических уравнений
- •Единственность решений
- •Метод последовательных приближений
- •Системы уравнений
- •Скоррелированные блуждания
- •Системы стохастических уравнений
- •Стохастический осциллятор
- •Линейные многомерные модели
- •Многомерие помогает одномерию
- •Как решать стохастические задачи?
- •Стохастическая природа
- •Теория броуновского движения
- •Стохастический осциллятор
- •Дрожание земной оси
- •Электронный шум
- •Хищники и их жертвы
- •Стохастическое общество
- •Финансовые рынки
- •Эмпирические закономерности
- •Диверсификация
- •Портфель на всю жизнь
- •Опционы
- •Кривая доходности
- •Компьютерное моделирование
- •Статистики
- •Случайные числа
- •Моделирование стохастических процессов
- •Ошибки вычислений и ускорение сходимости
- •Вычисление средних
- •R: Стохастический справочник
- •Основные соотношения теории
- •Системы уравнений с одинаковым шумом
- •M: Математические приложения
- •H: Помощь
- •C: Примечания
- •Рекомендуемая литература
244 |
Глава 9. |
9.3Случайные числа
При проведении численных моделирований нам необходимо уметь ге-
нерить случайные числа. Современные компьютеры в большинстве сво¼м являются цифровыми, а не аналоговыми устройствами. Поэтому возможна только генерация псевдослучайных чисел, являющихся результатом работы некоторого детерминированного алгоритма. В основе получения случайных чисел с различными распределениями лежит генератор равномерно распредел¼нных, обычно целых чисел из диапазона [0...max]. Понятно, что, если необходимо вещественное равномерно распредел¼нное число в диапазоне [0...1], то целое случайное число делится на максимальное значение max.
В C++ есть встроенная функция rand(), возвращающая целое, равномерно распредел¼нное случайное число в диапазоне от нуля до константы RAND_MAX, которая в большинстве компиляторов равна 32767.
Если в программе при помощи rand() генерить случайные числа, то при е¼ повторном запуске получится в точности такая же последовательность чисел. Чтобы е¼ изменить, можно в начале вызвать функцию srand(seed), где целое число seed зада¼т номер последовательности слу- чайных чисел. Если передать в эту функцию количество секунд, прошедших с 1970 г., возвращаемое функцией time(0), то последовательности будут получаться каждый раз разные.
Ниже в примере происходит генерация десяти случайных вещественных чисел из диапазона 0 6 x 6 1:
# include |
< time .h > |
// для функции t i m e |
|||
# include < stdlib .h > |
// |
для функций rand , s r a n d |
|||
# include |
< stdio .h > |
// |
для функции |
p r i n t f |
|
typedef |
double |
Float ; |
// |
вещественный |
òèï |
void main ( void ) |
|
|
|
|
|
{ |
|
|
|
|
|
srand ( time (0) |
% RAND_MAX ); |
// |
" встряхиваем " генератор |
||
for ( int i =0; |
i <10; i ++) |
|
|
|
printf (" %12.8 f\n" , Float ( rand ()) / RAND_MAX );
}
Если закомментировать srand, то при последовательных вызовах программы ряд случайных чисел будет получаться один и тот же. Операция % обозначает остаток от деления. Например, 47%10 равно 7. С е¼ помощью мы обрезаем значение seed величиной RAND_MAX-1. Для преобразования целых чисел в вещественные используется операция Float(rand()). Напомню, что 1/2 равно 0, тогда как Float(1)/2 равно 0.5.
Компьютерное моделирование |
245 |
В основе алгоритма rand() лежит формула, подобная следующей:
xt+1 = (16807 xt) mod (231 1): |
(9.5) |
Начиная с некоторого целого числа x0 6= 0, задаваемого при помощи функции srand(), при каждом вызове rand() происходит вычисление нового псевдослучайного числа на основе предыдущего. Магические константы, использующиеся в этом алгоритме, являются простыми числами, и качество генератора существенно от них зависит. Реализация этого алгоритма имеет вид:
unsigned int |
gRndSeed = |
1; |
// |
п о с л е д н е е случайное число |
inline void |
SRnd ( unsigned |
int |
seed ) // задание gRndSeed |
|
{ |
|
|
|
|
gRndSeed |
= ( seed ==0)? |
1: seed ; |
}
//
inline unsigned int |
RndI () // равномерное сл . число |
{ |
|
return gRndSeed |
= (16807 L* gRndSeed ) % 2147483647 L; |
} |
|
От хорошего генератора мы ожидаем выполнения четырех требований:
1)случайные числа равномерно распределены;
2)они независимы;
3)их последовательность имеет большой период повторения;
4)они различны.
Последние два требования не являются эквивалентными. Например плохой генератор может, имея формально большое RAND_MAX, выдавать только, скажем 1000 различных чисел, переставляя их самым разнообразным образом. При этом период повторения будет велик, но при построении вещественного генератора из всего континуума мы будем получать только 1000 случайных значений. Как мы увидим ниже, эта проблема присуща как библиотечной функции rand(), так и алгоритму (9.5).
Поэтому мы будем использовать генератор случайных чисел, основанный на эффекте переполнения 32-разрядных целых чисел. Соответствующая функция имеет вид:
inline unsigned int |
RndU () |
// равномерное сл . число |
{ |
|
|
return gRndSeed |
= gRndSeed *1664525 L + 1013904223 L; |
|
} |
|
|
Именно она будет в дальнейшем использоваться для генерации равномерно распредел¼нных случайных чисел.
246 |
Глава 9. |
Провед¼м сравнение статистических свойств тр¼х генераторов rand(),
RndI() и RndU(). Будем сразу работать с вещественными случайными числами в диапазоне [0...1]. Так, алгоритм переполнения RndU() имеет следующую вещественную реализацию:
inline Float Rnd () |
// равномерное сл . число [ 0 . . 1 ] |
{ |
|
return Float ( RndU ()) / |
Float (0 xffffffff ); |
} |
|
Константа 0xffffffff является записью максимального беззнакового целого 232 1 в шестнадцатеричной системе исчисления. Если целые
числа генерятся при помощи rand(), делить необходимо на RAND_MAX, à äëÿ RndI() íà 2147483647 = 231 1.
В качестве статистик, характеризующих равномерность случайных
чисел, мы будем использовать среднее, точное значение которого равноx = 1=2, квадрат среднего x2 = 1=3, минимальное xmin = 0 и макси-
мальное значение xmin = 1. Кроме этого, будем вычислять среднеквадратичное отклонение вероятности попадания в один из m одинаковых интервалов, на которые разобь¼м диапазон [0..1]:
|
n m |
1 |
|
|
2 |
||
zm = m k=1 |
pk m |
; |
|||||
|
|
X |
|
|
|
|
|
где n число сгенер¼нных случайных чисел, pk = nk=n, à nk-количество чисел, попавших в k-й интервал. Мы умножили эту статистику на n, так как подобное отклонение убывает, как 1=n, поэтому при любых n она будет примерно постоянной. Зная zm, можно по 2-критерию оценить достоверность гипотезы о равномерном распределении, однако мы будем использовать zm только для сравнения различных генераторов.
Характеристиками независимости будут служить ковариационíые коэффициенты для моментов случайных чисел, умноженные на pn:
cjk = pn xjt xkt 1 xjt xkt 1
В случае независимых величин чисел они должны быть равны нулю.
Для вычисления количества различных случайных чисел, выдаваемых генератором, воспользуемся библиотечным классом map который хранит пары ключ-значение . В нашем случае ключом будет служить случайное число, а значением количество его появлений. Класс map сравнительно быстро ищет и вставляет новые числа, и его размер size сообщает о количестве различных ключей (случайных чисел).
Компьютерное моделирование |
247 |
# include |
" stat . cpp " |
|
// |
статистическая |
библиотека |
||
# include |
<map .h > |
|
// |
map [ ] |
|
|
|
map < Float , unsigned int > |
lst ; |
|
|
|
|
||
const int |
n |
= 1000000; |
|
// |
величина выборки |
||
const int nex = 1000; |
|
// число экспериментов |
|||||
const int |
m |
= 100; |
|
// |
число |
интервалов |
|
void main ( void ) |
|
|
|
|
|
||
{ |
|
|
|
|
|
|
|
Float a[ nex ], s[ nex ], |
z[ nex ]; |
|
|
|
|||
Float c11 [ nex ], c22 [ nex ], c12 [ nex ], |
c21 [ nex ]; |
|
|||||
Float min [ nex ], max [ nex ], |
size [ nex ]; |
|
|
||||
Float p[m ]; |
|
|
|
|
|
|
|
for ( int ex =0; ex < nex ; |
ex ++){ // |
проводим 1000 |
экспериментов |
||||
srand ( ex +1); |
|
// |
" встряхиваем " |
генераторы |
|||
SRnd ( ex +1); |
|
|
|
|
|
||
a[ ex ]= s[ ex ]=0; |
|
// |
с р е д н е е и квадарат с р е д н е г о |
||||
c11 [ ex ]=0; |
|
// ковариации |
|
||||
c12 [ ex ]= c21 [ ex ]= c22 [ ex ]=0; |
|
|
|
||||
max [ ex ]=0; min [ ex ]=1; |
// |
максимальное и |
минимальное |
||||
Float rnd , rnd1 = Rnd (); |
// |
текущее и предыдущее число |
|||||
lst . clear (); |
|
// очищаем хэш таблицу |
|||||
for ( int k =0; k <m; k ++) p[k ]=0; |
|
|
|||||
for ( int |
i =0; i <n; |
i++ , rnd1 = rnd ){ |
|
|
|||
|
rnd = |
Rnd (); |
|
|
|
|
|
|
a[ ex ]+= rnd ; |
|
|
|
|
|
|
|
s[ ex ]+= rnd * rnd ; |
|
|
|
|
|
|
|
c11 [ ex ]+= rnd * rnd1 ; |
|
c12 [ ex ]+= rnd * Sqr ( rnd1 ); |
||||
|
c21 [ ex ]+= Sqr ( rnd )* rnd1 ; |
c22 [ ex ]+= Sqr ( rnd * rnd1 ); |
|||||
|
lst [ rnd ]++; |
|
|
|
|
|
|
|
int k= int ( rnd *m ); if (k >= m) k=m -1; if (k <0) k =0; |
||||||
|
p[k ]++; |
|
|
|
|
|
|
|
if ( max [ ex ]< rnd ) |
max [ ex ]= rnd ; |
|
|
|||
|
if ( min [ ex ]> rnd ) |
min [ ex ]= rnd ; |
|
|
|||
} |
|
|
|
|
|
|
|
a[ ex ]/= n; |
s[ ex ]/= n; |
|
|
|
|
|
c11 [ ex ]/= n; c11 [ ex ]=( c11 [ ex ]-a[ ex ]* a[ ex ])* sqrt ( Float (n )); c12 [ ex ]/= n; c12 [ ex ]=( c12 [ ex ]-a[ ex ]* s[ ex ])* sqrt ( Float (n )); c21 [ ex ]/= n; c21 [ ex ]=( c21 [ ex ]-a[ ex ]* s[ ex ])* sqrt ( Float (n )); c22 [ ex ]/= n; c22 [ ex ]=( c22 [ ex ]-s[ ex ]* s[ ex ])* sqrt ( Float (n ));
z[ ex ]=0; |
|
for ( int |
k =0; k <m; k ++) { z[ ex ]+= Sqr (( p[k ]/ n ) -1.0/ m ); } |
z[ ex ]*= Float (n )/ m; |
|
size [ ex ] |
= Float ( lst . size ())/ n; |
printf ("%d\n" , ex ); |
|
} |
|
printf ("%g\t%g\n" , Aver (a , nex ), Sigma (a , nex )/ sqrt ( nex )); // . . .
}
248 |
Глава 9. |
После получения на основании nex экспериментов массивов со ста-
тистиками при помощи функций Aver(), Sigma(), находящихся в файле stat.cpp, вычисляются среднее значение статистики и одна стандартная ошибка: Sigma(a, nex)/sqrt(nex). Ниже в таблицах стандартная ошибка используется для вывода только значащих цифр в статистике и указы-
вается в виде индекса. Так, например, 0.025 3 означает, что в последнем знаке может быть ошибка, т.е. 0:025 0:003. Привед¼м результаты рабо-
ты программы:
|
11 |
12 |
21 |
22 |
size=n |
RndU |
0.0033 |
0.0043 |
0.0033 |
0.0033 |
1.000000 |
rand |
0.0013 |
0.0013 |
0.0023 |
0.0013 |
0.032767 |
RndI |
0.252 |
0.273 |
0.273 |
0.283 |
0.049 |
Корреляционные свойства RndU и rand примерно одинаковы, а для RndIзаметно хуже. Генератор RndU существенно лучше по количеству различных чисел, которые он генерит. Так, на выборке из n = 1000000
он выда¼т все числа различными, тогда как для rand различных чисел только RAND_MAX=32767. Аналогичная проблема и у RndI, к тому же их количество существенно зависит от начального значения gRndSeed.
Если нам необходимо равномерно заполнить некоторый отрезок точка-
ми (например, при Монте-Карло вычислении интеграла), то, используя rand или RndI, мы получим около 105 точек, и проведение большего чис-
ла экспериментов становится бессмысленным. Статистики равномерности имеют следующие значения:
RndU |
|
1 |
|
|
1 |
0.009895 |
|
|
|
x |
|
x2 |
|
|
z100 |
min |
max |
|
0.49999 |
|
0.33332 |
|
|
0.000000 |
0.9999999 |
|
rand |
0.499991 |
0.333321 |
0.010095 |
0.000000 |
1.0000000 |
|||
RndI |
0.50041 |
|
0.33371 |
|
0.731 |
0.0000191 |
0.99998592 |
Средние значения у первых двух генераторов в пределах стандартной ошибки совпадают с точными значениями. Генератор RndI хуже как по этим показателям, так и по степени равномерности для m = 100 интер-
валов z100.
В дальнейшем для равномерно распредел¼нных чисел мы будем использовать RndU(). С помощью RndU(), а точнее, е¼ вещественного аналога Rnd(), можно получать и нормально распредел¼нные числа, которые, собственно, нам и необходимы для стохастических дифференциальных уравнений. Рассмотрим соответствующий алгоритм.
Компьютерное моделирование |
249 |
Для генерации нормально распредел¼нного случайного числа с нуле-
вым средним и единичной волатильностью служит метод Бокса-Миллера.
Перейд¼м от равномерно распредел¼нных на интервале [0::1] чисел x1 è p p
x2 к переменным: y1 = 2 ln x1 cos(2 x2), y2 = 2 ln x1 sin(2 x2). Количество чисел, попавших в некоторую область двухмерного пространства , определяется двойным интегралом. При замене переменных в
интеграле объ¼м умножается на якобиан:
Z |
Z |
|
|
|
|
Z |
|
|
|
|
|
|
|
|
|
|
|
@x1 |
@x1 |
|
p2 |
2 |
|
p2 |
2 |
|
|||||||
@y1 |
@y2 |
|
y1 |
|
|
y2 |
||||||||||
|
|
|
|
|
1 |
|
1 |
|
||||||||
|
|
@y1 |
@y2 |
|
|
|
e 2 |
|
|
|
e 2 dy1dy2: |
|||||
dx1dx2 |
= |
|
|
dy1dy2 = |
|
|
|
|
|
|||||||
|
|
|
@x2 |
@x2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Таким образом, переменные y1 è y2 являются независимыми (интегралы расщепляются) и имеют нормальные распределения.
Функции cos и sin вычисляются достаточно долго. Чтобы избежать
работы с ними, воспользуемся следующим при¼мом. Выберем область интегрирования в виде круга с единичным радиусом и центром в начале
координат. Будем случайным образом равномерно заполнять его точками. Пусть v1 è v2 координаты точки внутри круга v12 + v22 = r2 < 1.
В полярных координатах v1 = r cos , v2 = r sin . Òàê êàê v1 è v2 распределены равномерно, то и полярный угол будет равномерно распредел¼н в интервале [0::2 ]. Поэтому для вычисления значения cos и sin можно воспользоваться отношениями v1=r è v2=r:
Float RndG () |
|
|
|
|
// |
ã à ó ñ ñ î â î |
сл . число |
||
{ |
|
|
|
|
|
|
|
|
|
static |
int |
was |
= |
0; |
// |
была ли вычислена пара чисел |
|||
static |
Float |
r |
= |
0; |
// |
предыдущее случ . число из пары |
|||
if ( was ){ was =0; return |
r ;} // |
r из предыдущих вычислений |
|||||||
Float s , v1 , v2 ; |
|
|
|
|
|
||||
do { |
|
|
|
|
|
|
// |
жд¼м попадания в круг |
|
v1 |
= |
2* Rnd () -1; |
|
// |
точка в |
квадрате [ 1 . . 1 ] |
|||
v2 |
= |
2* Rnd () -1; |
|
|
|
|
|||
s= v1 * v1 + v2 * v2 ; |
|
// |
квадрат |
расстояния от центра |
|||||
} while (s >=1.0 |
|| |
s ==0); |
|
|
|||||
was |
= |
1; |
|
|
|
|
|
|
|
r |
= |
v1 |
* sqrt ( -2* log (s )/ s ); |
// п е р в о е число ( запомнили ) |
|||||
return |
v2 |
* sqrt ( -2* log (s )/ s ); |
// второе число ( вернули ) |
||||||
} |
|
|
|
|
|
|
|
|
|
Переменные was è r объявлены в функции RndG, как статические ( static). По завершении работы функции их значения будут сохранены, и при следующем вызове мы будем иметь значение одного из случайных гауссовых чисел, вычисленных на предыдущей итерации.