<<

стр. 4
(всего 5)

СОДЕРЖАНИЕ

>>

?1 9 1 3 2?
? ?
c. K = ?1 3 17 16 6 ? ,
?2 37 ?
6
? ?
?1 9 ?
??
m = ? 2 ?;
? 12 ?
??
4. Написать алгоритм моделирования однородного гауссовского процесса, заданного
начальным распределением
0 ? x ?1
? 2 x,
f ( x) = ?
?0, x < 0 или x > 1
и ковариационной функцией
?e ? x , x > 0
K ( x) = ? .
?0, x ? 0
5. Написать алгоритм моделирования одномерного стационарного гауссовского
марковского процесса с нулевым априорным средним и одномерной корреляционной



70
функцией вида
K (t , s ) = K (t ? s ) = ? 2 e ? b|t ?s| , b > 0.
6. Написать моделирующий алгоритм двумерного стационарного гауссовского
марковского процесса с нулевым априорным средним и матричной корреляционной
функцией вида
? 2 1?
K (t , s ) = K (t ? s ) = e ?b|t ?s| ?
? 1 1? , b > 0.
?
? ?
7. Написать моделирующий алгоритм одномерного нестационарного гауссовского
марковского процесса с нулевым априорным средним и корреляционной функцией
вида
? 3s
?t2 , t > s
K (t , s ) = K (t ? s ) = ? , t, s > 0.
3t
? 2, t?s
?s
Рассматриваемый процесс имеет место на выходе RC-фильтра при подаче на его вход
белого шума, если сопротивление R линейно изменяется во времени по закону
R = R(t) = t / 2.
R

C
Рис. RC-фильтр.

6.5. Контрольные вопросы и упражнения
1. Аппроксимация нестационарных случайных процессов.
2. Моделирование n-мерной случайной точки с независимыми координатами.
3. Моделирование n-мерной непрерывной случайной точки с зависимыми координатами.
4. Моделирование многомерного изотропного вектора.
5. Метод исключения.
6. Моделирование многомерного марковского гауссовского процесса.
7. Моделирование одномерного стационарного гауссовского процесса с дробно-
рациональной спектральной плотностью.
8. Моделирование процесса скользящего суммирования.
9. Моделирование нормального случайного процесса с использованием канонического
разложения.
10. Построить алгоритм моделирования псевдослучайного вектора ?, равномерно
распределённого на областях
a. симплекс
? ?
n
? = ? x ? R : ? xi = 1, x1 ,..., x n ? 0? ;
n

? ?
i =1

b. сфера
? ?
n
? = ? x ? R : ? xi2 = 1? ;
n

? ?
i =1

c. эллипсоид
? ?
n n
? = ? x ? R : ?? cij xi x j ? Q ? , где Cij – элементы положительно
n

? ?
i =1 j =1

полуопределённой матрицы; Q – заданная константа.

71
Глава 7

Общая схема метода Монте-Карло. Роль закона больших чисел и
предельных теорем в теории статистического моделирования.
Статистическая теория оценивания
Возможность моделирования случайных величин и процессов может быть использована
для моделирования реальных явлений и ситуаций. При этом наблюдение небольшого числа
реализаций случайной величины вряд ли принесет пользу, но наблюдение за большим их
числом позволяет сделать правильные выводы об их средних характеристиках. Такой подход
лежит в основе метода Монте-Карло, который использует различные предельные
соотношения теории вероятностей:
• закон больших чисел;
• предельные теоремы.

Пример 7.1.
Предложить схему вычисления числа ?=3,1415...
?=С/2L, где С - длина окружности, L - ее радиус. Рассмотрим окружность единичного
радиуса.
1


-1 1



-1

Рис. 7.1.

Генерируем две случайные величины ?1, ?2? Rav(0,1). Тогда координаты случайной
точки (?,?) вычисляются по формулам:
?=2?1-1
?=2?2-1
p=S?/S - вероятность попадания точки в круг.
S =4; S?=?r2=?; p=?/4
Таким образом ?=4р.
Пусть N - общее число точек, n - число точек попавших в круг. Тогда ˜ = n / N -
p
частотная оценка вероятности. Теперь оценку ? можно найти по формуле:
˜ = 4 ˜ = 1 4 ? , где ? = ?1, если точка в круге, то есть ? + ? ? 1;
? 2 2
N

?i
? p (7.1)
?
i
?0, если точка вне круга, то есть ? 2 + ? 2 ? 1.
N i =1 ?

Алгоритм:
1. генерируем пары точек (?1i,?2i) ? Rav(0,1);
2. проверяем неравенство (2?1-1)2+(2?2-1)2 ?1;
˜
3. вычисляем ? по формуле (7.1).

Уже при N=100 можно достичь хорошей точности оценки.


72
Теорема 7.1 (Теорема Колмогорова).
Для того, что бы среднее арифметическое независимых реализаций случайной величины ?
сходилось с вероятностью равной единице к ее математическому ожиданию, необходимо и
достаточно, чтобы это математическое ожидание существовало.

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

Пусть ? - произвольная случайная величина, М[?]=m.
1n
_ _
? n = ? ? i , ? n ??> m , при n > ? .
p =1
?1 , ?2 … ?n ,
n i =1

Рассмотрим примеры, где можно использовать вышеприведенные вычислительные схемы.

Пример 7.2.
Пусть требуется выполнить расчет надежности изделия, состоящего из значительного
числа элементов, для которых в результате эксперимента удалось построить распределение
времени безотказной работы.
Если это время распределено по показательному закону для всех элементов, то задачу
удастся решить аналитически. Если же имеют место отказы в результате износа элементов
распределенные по нормальному закону или усеченному нормальному закону, то
аналитическое решение невозможно.
Пусть изделие работоспособно, если работоспособны элементы, составляющие одну из
цепочек 1-7.

2 5


1 3 7
начало конец


4 6

Рис. 7.2.

Пусть ?i –случайная величина, равная времени безотказной работы i-ого элемента. Тогда
среднее время безотказной работы можно представить как математическое ожидание ?, где ?
=min{?1;max[min(?4,?6);min(max(?2,?3),?5)];?7}.

1n
Итак, закон больших чисел обеспечивает сходимость ? ? i к М[?].
n i =1

Пример 7.3.
Численное решение задачи Сильвестра.
Пусть задана выпуклая область D на плоскости и функция совместной плотности
распределения G( x1, y1; x2, y2; x3, y3; x4, y4 ) четырех точек Аi=(xi, yi)?D. Найти вероятность
того, что случайные точки образует выпуклый четырехугольник.




73
Рис. 7.3.
1, если полученный четырехугольник выпуклый;
Пусть ?=
0, в противном случае.

?1 n ?
Очевидно, что M ? ? ? i ? = p - искомая вероятность.
? n i =1 ?
Для зависимых случайных величин также имеют место законы больших чисел.
Пусть ?1 , ?2 , …, последовательность случайных состояний стационарной цепи Маркова с
переходной матрицей P.

Теорема 7.2 (усиленный закон больших чисел для цепи Маркова).
Если pT=|| p1, p2, …, pn || - единственное стационарное распределение стохастической
матрицы P, то для любого начального распределения p0 и произвольной функции состояний
f(?) имеет место равенство:
? ?
N n
?1
? ?
? f (? i ) = ? p j f (? j )? = 1,
P ?lim N
k
k =1 j =1
? ?
? ?
N >?

где i1>i2>… - последовательность переходов цепи Маркова.
Теорема 7.3 (закон больших чисел в ослабленной форме).
Для того, чтобы для последовательности ?1 , ?2 …- случайных величин таких, что
существует и конечно M[?i] ( при произвольном характере их зависимости) ? ?>0
выполнялось соотношение :
?1N ?
1N
P ? ? ? i ? ? M [? i ] < ? ? = 1
lim ? N i=1 N i=1 ?
N >?

необходимо и достаточно, чтобы выполнялось равенство:
? ?
N

? (?i ? M [?i ]) ?
2
?
? ?
lim ?=0
M ? i =1 N
? N + ? (?i ? M [? i ]) ?
2
N >? 2
? ?
? ?
i =1

Итак, законы больших чисел, хотя и могут служить принципиальным обоснованием
методов статистического моделирования, но они не позволяют судить об их погрешности.
Поведение погрешности можно изучить, используя предельные теоремы теории
вероятности.
Будем различать 4 основные задачи:

М[?2] –известные
?i – независимые;
А: конечные;
М[?2] –известные
?i – независимые;
Б: бесконечные;
М[?2] –известные
?i – зависимые;
В: конечные;
М[?2] –известные
?i – зависимые;
Г: бесконечные.


74
М[?] - неизвестно. Перед нами стоит задача вычисления математического ожидания
случайной величины ?.
В случае А для изучения поведения погрешности используется классическая предельная
теорема.
Пусть ?1 , ?2 , …, - независимые случайные величины, существуют и конечны: mk=M[?k] ,
n
? к=D[?k]; ? к= ? ? 2 k .
2 2

k =1


Теорема 7.4.
Если последовательность ?1 , ?2 , …, при любом постоянном ? > 0 удовлетворяет условию
Линдеберга:
1n
? ? ( x ? mi ) ?Fi ( x ) = 0,
2
lim ? 2 n i =1 (7.2)
N >? x ? mi < ?

где Fi(x) - функция распределения ?i ,то при n>? равномерно относительно x, получаем:

?1 ? x
n
1
? ?
2
e ? t / 2 ?t .
(?i ? mi ) < x ? >
P? (7.3)
? ?n 2? ??
?
i =1


Эта теорема сформулирована без предположения об одинаковой распределенности ?i.
При наличии его, условие (7.2) вытекает из конечности второго момента ?. Выражение (7.3)
в случае одинаково распределенных случайных величин можно записать в виде:
?1 n D(? ) ? x
? ? 1
P ? ? (? i ? mi ) < x ??e ?t .
?t / 2
?>
2? ?
? n i =1 n?
? ?

Считая n достаточно большим и оценку для D[?] достаточно точной, можно с заданной
вероятностью получить неравенства:
x ? D[? ] x ? D[? ]
1n 1n 2
M [? ] ? ? ? i ? M [? ] ? ? ? i ?
2

n i =1 n i =1
n n

x?
2
?
2
e ?t
P? = ?t ,
/2
если конечен четвертый момент.
? 0

Распределение среднего арифметического для задачи Б подробно изучено.
В случае существования r-го центрального момента порядок сходимости среднего
? 1? r ?
арифметического к M[?] оказывается равным O ? n r ? .
? ?
? ?
В случае В наиболее важна ситуация, когда случайные величины связаны в цепь Маркова.
Их можно рассматривать как реализации случайного процесса ?(t) в некоторые дискретные
моменты времени i ??t , i=0,1,…

Общий метод оценки математических ожиданий

Пусть ? - произвольная случайная величина.
М[?]=a – существует тогда, когда существует M[?2].
1
_ _
? N = ? ? i ; ? N ?P a , N>?
?1 , ?2 , … ?N; ?>
ni

75
т.к. последовательность одинаково распределенных независимых случайных величин, у
которых существует математическое ожидание, подчиняются закону больших чисел
(теорема Хинчина)
_
? N ? a , при больших N.

Определим погрешность метода.
Пусть ? имеет дисперсию D?=M(?2)-(M?)2.
Из курса теории вероятности известно, что последовательность одинаково
распределенных независимых случайных величин с конечными дисперсиями подчиняются
центральной придельной теореме, то есть ? х1 < x2
x2
? ?
N
1
P ? x 1 < (1 / ND ? ) ? (? i ? a ) < x 2 ? = ?
2
e ?t / 2 ?t
lim ? 2 ? x1
?
N>? i =1

Пусть x2=-x1=x , тогда:
?1 D? ?
x
2
(?i ? a) < x
lim? N ? ?
2
e ?t ?t
= ?( x) , ? ( x ) =
?
2? 0
N?
N >? ?

?_ ?
Следовательно, при больших значениях N: P ? ? N ? a < x D? / N ? ? Ф( x) .
? ?
Зафиксируем x? такой, чтобы Ф(x?)=?.
?_ ?
Тогда: P ? ? N ? a < x D? / N ? ? ? .
? ?

Эмпирическая оценка дисперсии:
2
?1 ?
1
1
? ??
?? ?i ? ? ?i ? .
2 2
2
М(? )? ; D??
i
N ?N
N

Для небольших N:
(? ? i )2 .
1 1
?
D? ? ?i 2 ?
N ?1 N ( N ? 1)

При использовании метода Монте-Карло моделируются случайные величины с
известными законами распределения, и из этих величин по заданным формулам
(алгоритмам) вычисляются значения более сложных случайных величин.
Закон распределения их не может быть найден аналитически, поэтому используются
средства математической статистики.

Различают два подхода: параметрический и непараметрический.

Параметрический подход предполагает, что искомое распределение известно с точностью до
значений конечного числа параметров.
_ _
? =|| ?1 , … ?m || , ??? , ? = ? (? 1 ,..., ? n ) - статистика.
Т



Свойства:
• несмещенность (асимптотическая несмещенность)
?_ ? ?_ ?
M ?? (?1 ,..., ? n )? = ? , M ?? (?1 ,..., ? n )? = ? ;
lim ?
? ? ?
n >?

• состоятельность;
76
? ?
_
P ? lim ? (? 1 ,..., ? n )= ? ? =1;
? n> ? ?
• эффективность (асимптотическая эффективность)
?_ ? ?_ ?
D ?? (?1 ,..., ? n )? = J , lim D ?? (?1 ,..., ? n )? = J ?1 .
?1

? ? ? ?
n >?

Часто дисперсия оценки является показателем ее качества. В общем случае могут быть
введены другие критерии:
а) точечное оценивание.
_
? (? 1 ,..., ? n ), ? (? 1 ,..., ? n ).
б) интервальное оценивание.
?
_
Если при некотором ? , 0 < ? < 1 ? (? 1 ,..., ? n ), ? (? 1 ,..., ? n ). могут быть выбраны так,
?

? _?
? ? ?? ,? ?
_
что , то интервал называют доверительным
P ?? <? <? ?
?=
?? ?
?? ?
интегралом для ?.
1. Метод наименьших квадратов (МНК),
2. Метод максимального правдоподобия (МП).
Пример 7.4.
_ _
Пусть ?i ? N( m,? ). Найдем m , ?
2 2

? ( x i ? m) 2 ? ?1 ?
n n n
1
?
L = ? f ( xi ,? ) = (2? ) ?
?
?n / 2 ?n
exp?? = exp ?? ( x i ? m) 2 ?
?
2? 2 ? (2? ) n / 2 ? n ? 2?
2
?
? i =1
i =1 i =1

? ?L 1n
_

? ?m = 0, m = n ? xi
? i =1
? _2
? ?L = 0, ? = 1
n _

? ( x i ? m) 2
? ?? n ? 1 i =1
2
?
1n
_
m = ? ? i - выборочное среднее
n i =1
_2
1n _

? ( xi ? m) 2 - выборочная дисперсия
?=
n ? 1 i =1
_
Можно показать, что m ? N( m; ?2/2 )
Изложенные выше методы оценивания среднего и дисперсии широко используются, когда
известно, что в результате моделирования получены независимые реализации случайной
величины ?.
Если ?1 ,…, ?n, полученные в результате моделирования, зависимы, то приходится
применять сложный аппарат. Последовательность {?i} рассматривают как значения
стационарного случайного процесса ?(t) в равноотстоящие моменты времени ?=0,1,.. Т .
(Ермаков С. , Михайлов Г. , 1926)
Непараметрический подход используется, когда параметрический вид плотности неизвестен,
т.е. когда априорная информация о характере результатов либо отсутствует, либо имеет
общий характер.
Например: f(x) – гладкая функция и т.д.
Характерная особенность непараметрических оценок – наличие систематической ошибки.
Методы оценивания плотности распределения с помощью ортогональных разложений для
обработки данных, полученных в результате моделирования были предложены А. Фроловым
и Н. Ченуевым.

77
Глава 8

Методы Монте-Карло

Содержание:
8.1. Простейший метод Монте-Карло....................................................................................78
8.2. Методы, использующие аналитические или классические численные процедуры
интегрирования................................................................................................................... 83
8.3. Метод существенной выборки........................................................................................ 84
8.4. Расслоенная выборка........................................................................................................ 86
8.5. Симметризация подынтегральной функции.................................................................. 88
8.6. Использование смещённых оценок................................................................................ 90

8.1. Простейший метод Монте-Карло для вычисления интеграла
Произвольная область G может быть:
• ограниченная или неограниченная;
• связанная или несвязанная (односвязанная или многосвязанная);
• выпуклая или невыпуклая.

xOy – плоскость
P=(x,y) – точка плоскости
dP=dxdy – элемент площади.

J= ? f ( p) P( p )?p , (8.1)
G


? P( p)?p = 1 .
где P(p) – заданная плотность вероятностей, определенная в G так, что
G

Следовательно, любой интеграл можно считать интегралом вида (8.1).

Площадь G ? SG , P1(p) =1/SG ? плотность распределения случайной точки равномерно
распределено в G.
Если f1(p)=SG?f(p) , то очевидно ? f ( p)?p = ? f1 ( p ) ? P1 ( p)?p .
G G1

Чтобы построить метод Монте-Карло для расчета (1) рассмотрим случайную точку Q с
плотностью P(p) и введем скалярную случайную величину Z=f(Q).
MZ = ? f ( p) P( p)?p = J
G



1. Рассмотрим задачу о приближенном вычислении интегралов.
J = ? f ( x) p( x)dx = ? f ( x)dF ( x) , ? p ( x) = 1 (8.2)
? ? ?

?1
n
J = ? f (? i )
n i =1
? i – независимые реализации случайной величины, имеющей функцию распределения
где
F (x ) .
Такой способ оценивания обладает следующими особенностями и преимуществами по
сравнению с классическими методами интегрирования:

78
1) естественным образом строится алгоритм интегрирования, если имеется алгоритм
моделирования случайной величины, распределенной по F ( x ) .
2) порядок убывания остатка

1n
? f ( x)dF ( x) ? ? f (? i )
n i =1
не зависит от размерности s пространства ?
3) последовательный характер оценивания интеграла; при увеличении размерности s
используются ранее вычисленные значения.
4) простая процедура оценки погрешности параллельно с основными вычислениями; так,
если сходимость распределения среднего арифметического к нормальному закону имеет
порядок:

?1
2),
O(n
то вычисляя среднее арифметическое квадратов значений функций в ? i мы можем оценить
дисперсию среднего арифметического и построить доверительный интервал,
соответствующий заданному уровню доверия.


? f ( x)dx по ограниченной области ? можно считать
Вообще любой интеграл J =
?
1
интегралом (1). Пусть площадь ? есть S ? . Причем p ( x) = . Тогда имеем:
S?
J = ? f ( x)dx = ? f1 ( x) p ( x)dx , где f1 ( x) = S ? f ( x) .
? ?



2. Рассмотрим задачу вычисления интеграла:


J = ? f ( x, y ) p ( x, y )dxdy
?


0 ? f ( x, y ) ? C . В трехмерном пространстве рассмотрим
Пусть в области ?
˜
цилиндрическую область: ? = ? ? (0; C )


й




Рис. 8.1. Область и поверхность

В этой области случайная точка имеет плотность:
(? ,? ) ? p( x, y ) Q (? ,? , µ )
˜
79
µ ? p( z )
1
˜
Р ( x, y.z ) = p ( z ) p ( x, y ) = p ( x, y )
C
˜
Определим Q следующим образом:
(? ,? ) ? p( x, y ) ˜
? Q (? ,? , µ ) .
µ
? p( z )
˜˜ ˜ ˜
Возьмем n реализаций: Q1 , Q2 ,…, Qn . Пусть v - количество точек Q , которые
оказались ниже поверхности z = f ( x, y ) . Тогда:
P(? = m) = Cn p m (1 ? p) n ? m , m = 0,1,..., n
m
˜
где P – вероятность того, что точка Q окажется ниже поверхности z = f ( x, y ) . Найдем эту
вероятность
f ( x, y )
˜ ( x, y, z )dz = 1 f ( x, y ) p ( x, y )dxdy = 1 J .
p = P( µ < f (? ,? )) = ? dxdy ? ?
p
C? C
? 0
J nJ
Таким образом p = , но M [? ] = np = , поэтому
C C
? C? , т.к M [ J ] = M ? C ? ? = C M [?] = C nJ = J .
?
J= ?n? n
n ?? nC
Эту оценку можно записать в виде привычном для метода Монте-Карло.
?C , если µ < f (? ,? )
Пусть случайная величина z = ? и если точкам z1 , z 2 ,…, z n
0 , если µ ? f (? ,? )
?
? r = 1 ? zi .
n
˜˜ ˜
соответствуют точки Q1 , Q2 ,..., Qn , то имеет место: J
n i =1
Это геометрический метод Монте-Карло.

3. Оценим точность этих двух методов.

Очевидно, что для существования D [ x ] = M ? x 2 ? ? ( M [ x ]) необходимо и достаточно,
2
??
чтобы существовал второй момент M ? x 2 ? .
??
Пусть f ( x, y ) ? L2 , L2 - множество функций f ( x, y ) для которых сходится интеграл:

? f (x, y ) p( x, y )dxdy .
2

?
1 n
? ? D[ f (? i ,?i )] , где
Тогда: D[ J n ] =
2
n i =1
D[ f (? ,? )] = ? f 2 ( x, y ) p ( x, y )dxdy ? [ ? f ( x, y ) p ( x, y )dxdy ]2 и
? ?
1 n
?
D[ J r ] = 2 ? D[ zi ] , где D[ zi ] = M [ z 2 ] ? ( M [ z ]) 2 , но
n i =1


80
M [ z 2 ] = C 2 P( µ < f (? ,? )) = CJ и
D[ z ] = CJ ? J 2 = C ? f ( x, y ) p ( x, y )dxdy ? [ ? f ( x, y ) p ( x, y )dxdy ]2 .
? ?
Очевидно, что если 0 ? f ( x, y ) ? C , то
2
? f ( x, y ) p( x, y )dxdy ? C ? f ( x, y ) p( x, y )dxdy
? ?

? ?
Следовательно: D[ J п ] ? D[ J r ] .
Последнее неравенство показывает, что в каком-то отношении простейший метод Монте-
Карло лучше геометрического. С точки зрения дисперсии усредняемой величины
простейший метод Монте-Карло всегда точнее геометрического.
? D[? ]
1n
Ранее было показано, что M [? ] ? ? ? i = rn ? .
n i =1 n
2
? ?t
1 2
dt ? ? = 0,6745 ,
Пусть P = = 2
?e
2 50
D[? ]
тогда rn = 0,6745 , где rn - наиболее вероятная ошибка, (8.3)
n
2
? 0,6745 ?
n = D[? ]? ?.
? rn ?
?
Итак, несмотря на то, что вероятная ошибка оценки J n простейшего метода меньше
?
вероятной ошибки оценки J r геометрического, при одинаковых n , но время, затрачиваемое
на расчет первой из них, может оказаться гораздо больше времени, затрачиваемого на расчет
второй.
Условимся говорить, что метод Монте-Карло построен, если:

? = 1 ??i
n
1) есть J
n i =1
?i
2) указаны формулы моделирования

Определение: Трудоемкостью алгоритма Монте-Карло называют произведение t * D[? ] ,
? , D[? ] - дисперсия с. в.
где t – время расчета одного значения

Трудоемкость зависит от:
- вида интегрируемой функции
- ЭВМ
- транслятора

Пример 8.1.
1 5
Вычислить интеграл: J = ? 5 x dx = .
6
0


81
?1
n
Пусть p( x ) = 1 . Тогда: J = ? 5 ? i , где f (? ) = 5 ? i
n i =1
2
1
?5? 5
5
D[5 ? i ] = ? x 2 dx ? ? ? = .
?6? 252
0
Так как 0 ? 5 x ? 1 , то воспользуемся геометрическим методом.

?1
n
J = ? zi
n i =1


?1 , ? i1 < 5 ? 2i
?
zi = ?
?0 , ?1i ? 5 ? 2i
?

2
5 ?5? 5
, D[5 ? ] < D[ z ]
D[ z ] = ? ? ? =
6 ?6? 36

Оценим количество операций, затрачиваемых на расчет одного значения 5 ? , z .
? ? , вычисление
5 5
...
1) для вычисления - три операции для через
ln x, e x ? 70 операций = t
z - 6 операций для ? , проверка неравенства: (? 1i ) 5 ? ? 2i ? 10
2) для вычисления
операций = t

tD[5 ? ] ? 1,4
tD[ z ] ? 1.4 .

Однако простейший метод можно реализовать с помощью другого метода:
? = 1 ? max(? 1i , ? 2i , ? 3i , ? 4i , ? 5i ) ,
n
J
n i =1
тогда t = 20 операций и tD[max(...)] = 0.40 .

D [? ]
Из (8.3) видно, что вероятная ошибка пропорциональна . Скорость убывания этой
n
ошибки с ростом n невелика.
Важно научиться выбирать для расчета интегралов такие схемы, такие случайные
величины, для которых D [? ] по возможности мала.
Общие методы понижения дисперсии в целом дают значение дисперсии, меньшее чем у
простейшего метода Монте-Карло.

82
8.2. Методы, использующие аналитические или классические
численные процедуры интегрирования
Они предполагают, что часть задачи можно решить аналитически или численно.

а) Выделение главной части
Это общий принцип, относящийся ко всем методам Монте-Карло: если главную часть
задачи можно вычислить аналитически, то лучше считать методом Монте-Карло не всю
задачу, а только поправку-разницу между всей задачей и главной частью.
Пусть надо вычислить J = ? f ( x) p ( x)dx, f ( x) ? L2 и имеется функция g ( x) ? L2 “близкая”
X

? g ( x) p( x)dx = C
к f ( x) , такая, что значение интеграла известно.
X
n
1
Тогда J ? = C + ? [ f (?i ) ? g (?i ) ] ?i ? p( x) .
?
n i =1
И если ? = f ( x) ? g ( x) , то
? ??
2
?
1? ?
?
D [ J ?] = ? ? ? 2 ( x) p ( x)dx ? ? ? ?( x) p ( x)dx ? ? .
n ?X ?X ??
? ?
? [ f ( x) ? g ( x)]
2
p( x)dx ? ? , то очевидно, что
Если g ( x) настолько близка к f ( x) , что
X
1
?
D [ J ?] ? ?.
n
Общие правила выделения главной части указать трудно: в различных задачах это
делается по-разному.

б) Интегрирование по части переменных (понижение порядка интеграла)
Если аналитически взять интеграл по некоторым из переменных, а по остальным
переменным использовать метод Монте-Карло, то дисперсия уменьшится.
Пусть требуется вычислить J = ?? f ( x, y ) p ( x, y )dxdy .
X ?Y
Предположим, что аналитически можно произвести интегрирование в выражениях для
p1 ( x) и f ( x) , где
? p1 ( x) = p( x, y )dy
?
? Y
, т.е. предположим, что эти интегралы
?
? f ( x) = ? f ( x, y ) p( x, y ) [ p1 ( x)] dy = M ? f ( x, y )
?1
x?
? ?
?
? Y

можно вычислить аналитически.
˜
Тогда исходный интеграл запишется как J = ? f ( x ) p1 ( x )dx .
X
Теорема 8.1.
n
? ?? = 1 ? f (? ) где ? - распределена с
При оценке интеграла J с помощью оценки J i i
n i =1
плотностью p1 ( x) , дисперсия будет разве лишь меньше, чем при его оценке с помощью
n
? = 1 ? f (? , ? ) , где (? , ? ) распределена с плотностью p ( x, y ) .
оценки J i i i i
n i =1


83
Доказательство.
1? ?
1
?
1) D [ J n ] = 2 n D [ f (?, ?)] = ? ?? f 2 ( x, y ) p ( x, y )dxdy ? J 2 ?
n n ? X ?Y ?
1? ?
1
?
2) D [ J ??] = 2 n D[ f (?i )] = ? ? f 2 ( x) p1 ( x)dx ? J 2 ?
n n ?X ?
? ?
1? 2
3) D [ J n ] ? D [ J ??] = ? ? f ( x, y ) p( x, y ) ? f ( x) p1 ( x) ? dx =
? ? 2

n X ?Y ?
? ?
p1 ( x ) p ( y x )

? ??
2
?
? ?
? p ( x, y ) ? ?
1
= ? p1 ( x) ? ? f ( x, y ) p( y x)dy ? ? ? f ( x, y ) dy ? dx =
2
?Y ??
?Y
nX p1 ( x)
? ??
?
? ??
?
? ?
p( y x)


1 1
p1 ( x) ? M [ f 2 ( x, y ) x] ? ( M [ f ( x, y ) x]) ?dx = ? D[ f ( x, y ) x] p1 ( x) dx ? 0
2

n?
=
? ? n
X ?0
?0
p ( x, y )
Распределение y есть , где x - фиксировано.
? p( x, y)dy
Y
Итак интегрирование по части переменных может уменьшить дисперсию оценки, но
˜
нужно иметь в виду, что f ( x) может быть значительно более сложной функцией, чем f (x)
?
и следовательно количество работ по вычислению J ?? может существенно превосходить
?
количество работ по вычислению J n .

8.3. Метод существенной выборки (выборки по важности)
Рассмотрим общую схему метода.
J = ? f ( x) p ( x)dx , где ? p ( x)dx = 1 .
Пусть плотность g ( x) > 0 для тех x , для которых p ( x) > 0 . Тогда очевидно, что
p( x)
J = ? f ( x) g ( x)dx
g ( x)
? ??? = 1 ? f (? i ) p (? i )
n
? i ? g ( x)
J (8.4)
g (? i )
n i =1
? ?
J ??? = J ???( g ) образуют семейство несмещенных оценок интеграла J .
Оценки
Семейство связано с возможностью выбора g (x) . Можно поставить задачу о выборе
оптимальной плотности g = g 0 при которой дисперсия минимальна.

Теорема 8.2.
Минимальная дисперсия реализуется для плотности g 0 ( x) пропорциональной f ( x) p ( x)

{ }
1? 2
f ( x) p ( x)dx ? ? J 2
?
n ??
и выражается равенством D [ J ???( g 0 )] = (8.5)
?

Доказательство.


84
Если g 0 ˜ f ( x) p ( x) , тогда g 0 ( x) будет тогда плотностью, когда
f ( x) p ( x)
g 0 ( x) = . (8.6)
? f ( x) p( x)dx
? p (? i ) ?
n
1
?
?
И D [ J ???( g 0 )] = 2 D ? f (?i ) ? , но
g (? i ) ?
n ?
i =1

? p (?i ) ?
? = просто без индекса ? i = 1, n =
D ? f (?i )
g (?i ) ?
?
? p ( ?) ? p 2 ( x)
? = ? f ( x) g 2 ( x) g ( x)dx ? J =
= D ? f ( ?) 2 2

g (?) ?
?
f 2 ( x) p 2 ( x)
= подставим вместо g (x) выражение (8.6) = ? dx ? f ( x) p ( x)dx ? J 2 =
f ( x) p( x)
2
= ? ? f ( x) p ( x)dx ? ? J 2
? ?
?
Отсюда следует вывод: интеграл f ( x) p( x)dx должен быть конечным.
Оптимальность плотности g 0 ( x) следует из:
2 2
? ???( g )] ? D [ J ???( g )]? = f ( x) p ( x)dx ? J 2 ? ?
2
f ( x) p( x)dx ? + J 2
?
? g ( x) ?
?D [J
? 0? ? ?
? f ( ?) p ( ?) ?
=D ? ??0
g ( ?) ?
?
где ? распределена с g ( x) .
Из (8.6) следует, что для знакопостоянной функции оптимальная плотность требует для
своего построения знания точного интеграла J = ? f ( x) p ( x)dx . Вычисление его может быть
не менее сложной задачей, чем исходная.
Если выбирать g ( x) так, чтобы она грубо походила на f ( x) p ( x) , то можно надеяться
получить не слишком большое значение дисперсии.

Замечание. Метод существенной выборки позволяет обеспечить конечную дисперсию при
оценивании с помощью (8.4) интегралов от абсолютно интегрируемой с весом p( x) функции
f ( x) .

Замечание. Иногда возникает задача вычисления m интегралов J j = ? f j ( x) p( x)dx.


? ??? = 1 ? f (? ) p (?i ) , j = 1, m
n
Jj j i
g (?i )
n i =1
Пусть заданы некоторые вещественные константы a j , j = 1, m .

? ???] = min 1 ? ? a 2 f 2 ( x) p ( x) dx ? ? a 2 J 2 ? .
2
m m m
g ( x) ? min ? a D [ J j
n ? ? j =1
2
? j j?
j jj
g ( x)
g
?
j =1 j =1

Минимум достигается тогда, когда достигается минимум выражения:
p 2 ( x)
m
S = ? ? a j f j ( x)
22
dx .
g ( x)
j =1




85
1
? ? p ( ?)
m 2
Если обозначить ? = ? ? a 2 f j2 (?) ? , то S = ? ?2 g ( x)dx .
j
? g ( ?)
? j =1
1
? ?
m 2
D [?] = S ? J 2 , где J 2 = ? ? ? a 2 f j2 ( x) ? p( x)dx .
j
? j =1 ?
2
?m ?
1
? ?2
Так как D [?] = S ? J 2 ? 0 , то S ? ? ? ? ? a 2 f j2 ( x) ? p( x)dx ? . (8.7)
? j =1 j ?
? ?
? ?
? ?
2
?m ?
1
? ?2
Нижняя граница в (8.5) при S = ? ? ? ? a 2 f j2 ( x) ? p ( x)dx ? = J 2 ? min S ( g ) = J 2 ,
? j =1 j ? g
? ?
? ?
? ?
? ???] = 1 ? J 2 ? ? a 2 J 2 ? .
m m
а следовательно минимум выражения равен: min ? a j D [ J j 2
? j j?
n? ?
j =1 j =1
1
? ?
m 2
?? j a 2 f j2 ( x) ? p ( x)
? j =1 ?
Согласно теореме это будет при g ( x) = g 0 ( x) = .
1
?m 2 2 ? 2

? ? ? a j f j ( x) ? p( x)dx
? j =1 ?
Проверяется непосредственной подстановкой g 0 ( x) в S вместо g (x) .
Полученный результат важен, когда моделирование каждой из плотностей, близких к
f j ( x) p( x) - достаточно трудоемкая задача, и целесообразно ограничиться оптимальным
выбором одной плотности, “средней” между оптимальными плотностями.

8.4. Расслоенная выборка (или выборка по группам)
Рассматриваемый прием близок к методу существенной выборки.
Представим область, по которой производится интегрирование, в виде суммы конечного
числа непересекающихся подобластей Ai , так что
m

? f ( x) p( x)dx = ? ? f ( x) p( x)dx .
i =1 Ai
A

Каждый из интегралов может быть оценен с помощью простейшего метода Монте–Карло
p ( x)
? f ( x) p( x)dx = S ( Ai ) ? gi ( x) dx
S ( Ai )
Ai
ni
1
Этот интеграл оценивается с помощью суммы S ( Ai )? gi ( xki )
ni ki =1

? f ( x), x ? Ai
, S ( Ai ) = ? p( x)dx
где gi ( x) = ?
0 , x ? Ai
? Ai

xi с
xki ? (ki = 1,… , ni , i = 1,… , m) - реализация независимых случайных величин
p( x)
, x ? Ai .
плотностью вероятности
S ( Ai )



86
ni
m
1
Таким образом ??[ f ] = ? S ( Ai )? gi ( xki ) ? ? f ( x) p ( x)dx
i =1 ni ki =1

Очевидно, что M ??[ f ] = ? f ( x) p ( x)dx , а
A
2
1? ?
m m
1
D ? [ f ] = ? S ( Ai ) ? f 2 ( x) p( x)dx ? ? ? ? f ( x) p( x)dx ?
?
i =1 ni i =1 ni ? Ai ?
? ?
Ai



Пример 8.2.
n1 = p1 N ?
Пусть m = 2, p1 = S ( A1 ), p2 = S ( A2 ), n1 + n2 = N , ? - целые числа.
n2 = p2 N ?
1? 2 a 2 b2 ?
Тогда D ? = ? ? f ( x) p ( x)dx ? ? ? ,
?
N? p1 p2 ?
a = ? f ( x) p ( x)dx,
A1
(a + b) = ? f ( x) p( x)dx
где
? f ( x) p( x)dx,
b= A
A2

Для обычного метода Монте-Карло:

{ }
1
N?
D ?N [ f ] = f 2 ( x) p ( x)dx ? a 2 ? b 2 ? 2ab ? (a + b) 2
?

12 12
a + b ? ( a + b) 2
Легко показать, что (8.8)
p1 p2
Откуда D ? ? D ? для любой интегрируемой с квадратом относительно веса p (x)
?
функции f .
?1 ? ?1 ?
(8.8) равносильно неравенству: ? ? 1? a 2 + ? ? 1? b 2 ? 2ab ? 0 .
? p1 ? ? p2 ?
2
? ?
? a 1 ? 1 ? b 1 ? 1 ? ? 0 , т.к. p1 + p2 = 1 .
Следовательно, ? ?
p1 p21
? ?
Сравниваемые дисперсии равны в том случае, если aS ( A1 ) = bS ( A2 ) .
Итак – любое разбиение исходной области интегрирования при выборе n1 и n2
пропорциональном S ( A1 ) и S ( A2 ) ведет к уменьшению дисперсии, по сравнению с оценкой,
где выборка производится равновероятно из всей области.
Если n1 , n2 выбирать не пропорционально интегралам от p( x) , т.е. S ( A1 ) и S ( A2 ) , то
этот факт не имеет места.
Пример 8.3.
3 1
Пусть m = 2, n1 = n2 , S ( A1 ) = , S ( A2 ) =
4 4
? ??
2 2
? ??
1 ?3 ?
1
D ? = ? ? f 2 ( x) p ( x)dx + ? f 2 ( x) p( x)dx ? ? ? f ( x) p ( x)dx ? ? ? ? f ( x) p( x)dx ? ? .
?
?A ? ?A ?
n1 ? 4 A1 4 A2 ??
?1 ? ?2
? ?




87
Сравнивая полученное выражение с дисперсией обычного метода Монте-Карло с числом
1? 2 ??
2
?
? ?
? ? f ( x) p( x)dx ? ? ? f ( x) p( x)dx ? ? , убеждаемся, что есть
испытаний в 2n2 : D ? =
2n1 ? A ?A ??
? ?
функции, для которых такая расслоенная выборка дает худший результат, чем оценка при
помощи обычного метода Монте-Карло.
Для оптимального выбора n1 и n2 следует привлечь информацию о вторых моментах
оценок.
? f ( x), x ? A1 ? f ( x), x ? A2
Пусть g1 ( x) = ? g 2 ( x) = ?
0 , x ? A1 ?0 , x ? A2
?
gi = g (? i ), i = 1, 2, ? ? Rav( Ai )
2 2
1? ? 1? ?
?? ??
D ? [f ]= p ( x)dx ? Dg1 + p( x)dx ? Dg 2 ? min
?
p1 N ? A1 ? p2 N ? A2 ?
? ? ? ?
Так как p1 + p2 = 1 , то
Dg1 S ( A1 )
p1 =
S ( A1 ) Dg1 + S ( A2 ) Dg 2

[ ]
1 2
Dmin = S ( A1 ) Dg1 + S ( A2 ) Dg1
N
Является основой для построения многоэтапных методов Монте-Карло.

8.5. Симметризация функции
Различают простую и сложную симметризацию.
а) Простая симметризация
b b
1
J = ? f ( x)dx =(b ? a) ? f ( x) p( x)dx , где p( x) = .
b?a
a a

Простейший метод Монте-Карло приводит к оценке:
? = b ? a ? f (? ) = 1 ? z , z = (b ? a) f (? ) .
n n
J i i
i
n i =1 i
n i =1
1
( 1)
Рассмотрим симметризованную функцию: f ( x) = [ f ( x) + f (a + b ? x)]
2
?c = b ? a ? [ f (? i ) + f (a + b ? ? i )] = 1 ? zi(1) ,
n
n
J
2n i =1 n i =1
1
zi(1) = (b ? a) [ f (? i ) + f (a + b ? ? i )] - симметризованная оценка.
2
?
Можно показать, что M [ J c ] = J , а в каком соотношении находятся дисперсии?
b
M [ z ] = (b ? a ) ? f 2 ( x)dx .
2

a
b

? f ( x) p ( x ) dx )
2
( M [ z ] = (b ? a ) M [ f (?)] = (b ? a )
2 2 2 2
(8.9)
a
1
(b ? a )


88
(b ? a ) b 2
(1) 2 2
]= ? [ f ( x ) + 2 f ( x) f ( a + b ? x ) + f (a + b ? x)]dx =
M [z
4a
(b ? a) b 2 b
[ ? f ( x)dx + ? f ( x) f (a + b ? x)dx]. (8.10)
2 a a
Легко показать, что M [ z (1) 2 ] ? M [ z 2 ] , а следовательно D[ z (1)2 ] ? D[ z 2 ] .
(1 )
надо вычислить два значения f ( x ) .
Однако, для расчёта одного значения z
? ?
Поэтому трудоёмкость оценки J
c будет меньше трудоёмкости J только тогда, когда
D [ z (1 ) ] по крайней мере вдвое меньше, чем D[ z ] .
Для монотонных функций это всегда выполнено.

Теорема 8.3.
f ( x ) монотонна при a ? x ? b , то
Если кусочно-непрерывная функция
1
( 2 D [ z (1 ) ] ? D [ z ] ? 0 ) .
D[ z (1) ] ? D[ z ]
2
Доказательство.
b
D[ z ] = (b ? a) ? f 2 ( x)dx ? J 2 .
a
b b
2 D[ z ] = (b ? a) ? f ( x)dx + (b ? a) ? f ( x) f (a + b ? x)dx ? 2 J 2
2
(1)


a a
Утверждение теоремы равносильно неравенству
b
(b ? a ) ? f ( x) f (a + b ? x)dx ? J 2
a
Пусть для определённости f (x) не убывает и f (b) ? f (a) .
Введём вспомогательную функцию.
x
V ( x) = (b ? a) ? f (a + b ? t )dt ? ( x ? a) J
a
V (a ) = V (b) ? 0
V ?( x) = (b ? a) f (a + b ? x) ? J - монотонная (8.11)
V ?(a) > 0
, следовательно, V ( x) ? 0 при a ? x ? b .
?(b) < 0
V
b
Можно записать очевидное неравенство ? V ( x) f ?( x)dx ? 0 .
a
b
Проинтегрировав его по частям, получим неравенство ? f ( x)V ?( x)dx ? 0 , подставляя
a
сюда выражение (8.11) получаем, что требовалось доказать.
Случай не возрастания f (x) рассматривается точно также.


89
б) Сложная симметризация
Интервал (a; b) можно разбить на конечное число частей и на каждой из них использовать
простую симметризацию. Рассмотрим случай разбиения (a; b) на две равные части.
1
Пусть c = (a + b) . Тогда
2
b c b1c 1b
J = ? f ( y)dy = ? f ( y)dy + ? f ( y)dy = ? [ f ( y) + f (a + c ? y)]dy + ? [ f ( y) + f (c + b ? y)]dy
2a 2c
a a c
x = 2y ? a (a ; c) ? (a ; b)
В первом интервале сделаем замену переменной
x = 2y ?b (c ; b) ? (a ; b)
во втором интервале
b a+x 2a + b ? x 2b + a ? x
1
( 2)
Получаем J = ? f ( x)dx , где f (2) ( x) = [ f ( )+ f ( )+ f ( )]
4 2 2 2
a
1 n
?
J сл.с = ? zi( 2) , где zi( 2) = (b ? a ) f ( 2) (? i ) , ? i ? Rav(a; b) .
n i =1
Замечание. Различные методы симметризации, весьма наглядные и эффективные в
одномерном случае, становятся громоздкими и трудно оцениваемыми при переходе к
функциям многих переменных.

Замечание. Симметризация f ( x, y, z ) по всем переменным в единичном кубе содержит 8
слагаемых.
1
f (1) ( x, y, z ) = [ f ( x, y, z ) + f (1 ? x, y, z ) + f ( x,1 ? y, z ) + f ( x, y,1 ? z ) + f (1 ? x,1 ? y, z ) +
8
f (1 ? x, y,1 ? z ) + f ( x,1 ? y,1 ? z ) + f (1 ? x,1 ? y,1 ? z )]

<<

стр. 4
(всего 5)

СОДЕРЖАНИЕ

>>