Лекция 12. Метод статистических испытаний

1. Вычисление определённого интеграла

Пример 1.   Рассмотрим определенный интеграл

a

Δ =

1  ∫  0

g(x) dx,

где g(x) - ограниченная на отрезке [0,1] функция (|g(x)| ≤ c < +∞ для всех x О [0,1]). Рассмотрим также СВ X, равномерно распределенную на отрезке [0,1], для которой f_X(x) = 1, если x О [0,1] и f_X(x) = 0, если x П [0,1]. Тогда очевидно, что

M[g(X)] =

+∞  ∫ -∞

g(x)fX(x) dx =

1  ∫  0

g(x) dx

Δ =

a .

Таким образом, значение определенного интеграла можно найти как математическое ожидание СВ

~ X

Δ = g(X),

где СВ X имеет равномерное распределение на отрезке [0,1], X ~ R(0,1). Пусть теперь последовательность {X_n}, n = 1,2,... состоит из СВ, независимых и равномерно распределенных на отрезке [0,1]. Можно убедиться, что в этом случае СВ

~ Xn	Δ = g(Xn), n = 1,2,...,

будут также независимы, одинаково распределены и с ограниченной дисперсией. Поэтому по усиленному закону больших чисел (теорема Колмогорова) последовательность СВ

Yn

Δ =

1 n

n  ∑  k=1

~ Xk

сходится почти наверное к величине a = M[g(X)], т.е.

Yn

п.н. →

a.

Это значит, что, проведя большое количество испытаний, можно с высокой точностью вычислить значение интеграла a.

Определение 1.   Метод вычисления некоторой детерминированной величины a как среднего арифметического Y_n независимых СВ X_n с одинаковым распределением, подобранным таким образом, чтобы

Yn

п.н. →

a,

называется методом статистических испытаний (методом Монте-Карло) для вычисления a.

Замечание 1.   В некоторых учебниках метод Монте-Карло называют также методом статистического моделирования, что, на наш взгляд, не совсем корректно. Дело в том, что термин "моделирование'' используется также в теории математического моделирования для описания процесса создания математических моделей каких-либо явлений. Подчеркнем, что в методе Монте-Карло речь идет о статистической имитации испытаний (опытов), а не о процессе создания статистических моделей опытов. Оправдано говорить лишь о моделировании случайных величин, если используются, например, датчики случайных чисел.

Замечание 2.   Метод Монте-Карло имеет огромную область приложений. Наиболее трудной проблемой в его реализации является выбор необходимого числа испытаний n такого, чтобы можно было считать, что СВ Y_n достаточно "близка" к a. Ясно, что из-за вычислительных трудностей желательно выбирать величину n с возможно меньшим гарантирующим значением N. Для выбора N обычно используют центральную предельную теорему, считая, что распределение нормированной суммы Z_N СВ X_k, k = 1,N, является стандартным нормальным распределением N(0,1). Рассмотрим на примерах, как выбирается величина N.

Пример 2.   Выберем гарантирующее число статистических испытаний N при вычислении значения a определенного интеграла из примера 1. Рассмотрим дисперсию СВ

~ X

Δ = g(X),

предполагая, что g(x)

/ ≡

a и g(x) - непрерывна,

D[g(X)] =

+∞  ∫ -∞

[g(x)-a]2fX(x) dx =

1  ∫  0

[g(x)-a]2dx ≤

1  ∫  0

[|g(x)|+|a|]2dx .

По условиям примера 1 имеем |g(x)| ≤ c, тогда

|a| = |

1  ∫  0

g(x) dx| ≤ c ,   D[g(X)] ≤

1  ∫  0

[|g(x)|+c]2dx ≤ 4c2.

Учтем, что к последовательности СВ g(X_n), n = 1,2,..., применима центральная предельная теорема, в соответствии с которой СВ

Zn

Δ =

1 sn

n  ∑  k=1

[

g(Xk) - mk

]

, где sn2

Δ =

D

[

n  ∑  k=1

g(Xk)

]

, mk

Δ =

M[g(Xk)],

сходится по распределению к СВ U, т.е.

Zn

F →

U.

Убедимся в том, что в данном случае применима центральная предельная теорема, т.е. выполнено условие Ляпунова. С этой целью, учитывая, что СВ X_k, k = 1,n, имеют одно и тоже распределение R(0,1), найдем

mk

Δ =

M[g(Xk)] = M[g(X)] = a,

sn2

Δ =

D[

n  ∑  k=1

g(Xk)]

4)M[X]    =

n  ∑  k=1

D[g(Xk)] = nD[g(X)].

Тогда условие Ляпунова выполнено:

1  sn3

n  ∑  k=1

M[|g(X) - a|3] ≤

1 (nD[g(X)]))3/2

n  ∑  k=1

(2c)3 ≤

(2c/ε)3/2     √n

→ 0,

т.к. g(x)

/ ≡

a, g(x) - непрерывна и поэтому D[g(x)] ≥ ε > 0. Тогда

Zn = [

1 n

n  ∑  k=1

g(Xn) - a ]

√

n       D[g(X)]

= (Yn - a)

√

n       D[g(X)]

,

где

Yn

Δ =

1 n

n  ∑  k=1

g(Xk).

Далее, по усиленному закону больших чисел, СВ Y_n сходится почти наверное к постоянной

a

Δ =

M[g(X)],

т.е. при больших n СВ Y_n "близка" к a. Предположим, что нужно вычислить a с точностью Δ: |Y_n-a| ≤ Δ. Выясним степень достоверности этой оценки (Y_n является случайной величиной):

P{|Yn - a| ≤ Δ} = P

{

|Yn - a|

√

n       D[g(X)]

≤ Δ

√

n       D[g(X)]

}

.

Так как D[g(X)] ≤ 4c², то

P{|Yn - a| ≤ Δ} ≥ P

{

|Yn - a|

√

n       D[g(X)]

≤ Δ

√n 2c

}

= P

{

|Zn| ≤ Δ

√n 2c

}

.

В соответствии с центральной предельной теоремой,

Zn

F →

U, т.е.

P

{

|Zn| ≤ Δ

√n 2c

}

≈ 2Φ0[Δ

√n 2c

],

где

Φ0(x)

Δ =

1 √2π

x  ∫  0

exp[-

x2  2

] dx.

Выберем некоторую величину β (доверительная вероятность), которая равна требуемой надежности выполнения неравенства |Y_n - a| ≤ Δ, т.е. P{|Z_n| ≤ Δ√n / 2n} = β. На практике β обычно задают в пределах от 0.95 до 1. Тогда, решая уравнение 2Φ₀(Δ√N/(2c)) = β относительно N можно найти гарантирующее число испытаний N_β. Функция Ф₀(x) задается таблично, и поэтому для доверительной вероятности β можно найти такое значение x_β/2, что Ф₀(x_β/2) = β/2. Таким образом, получаем уравнение Δ√N / (2c)= x_β/2, из которого легко найти искомое значение

Nβ =

[(

2cxβ/2    Δ

)

2

]

+ 1,

где символ [.] означает целую часть дробного числа. Подчеркнем, что для вычисления N_β необходимо задать точность Δ оценки величины a и уровень доверительной вероятности β.

Замечание 3.   В данном примере рассмотрен случай вычисления

a

Δ =

M[g(X)]

для равномерно распределенной СВ X ~ R(0,1). Аналогично можно поступить при вычислении

a

Δ =

M[g(X)]

для произвольной СВ X, такой, что a < +∞. Тем самым дан ответ на вопрос (см. замечание Л10.Р2.З5), в каком смысле понимать точность приближения Y_n ≈ m_X.

2. Вычисление вероятности события

Пример 1.   Рассмотрим еще одно применение метода статистических испытаний для оценивания неизвестной вероятности

p

Δ =

P(A)

некоторого события A. В соответствии с теоремой Бернулли

M / n

п.н. →

p ,.

а в соответствии с теоремой Муавра - Лапласа

* Wn = Zn

Δ =

M-np √npq

F →

U, где U ~ N(0,1).

Воспользуемся этими теоремами для выбора гарантирующего числа испытаний N, при котором можно было бы сказать, что СВ M / N "близка" к оцениваемой вероятности p. Поступим аналогично тому, как выбиралось N_β в примере 2 из предыдущего раздела. Вначале зададим точность Δ оценки вероятности p: |M/n - p| ≤ Δ. Затем найдем надежность этой оценки, учитывая, что

* Wn = Zn ≈ U,

P

{

|

M  n

- p| ≤ Δ

}

= P

{

M-np √npq

≤ Δ

√

n   pq

}

≈ 2Φ0

(

Δ

√

n   pq

)

.

Задав доверительную вероятность β, получим трансцендентное уравнение: 2Φ₀(Δ√N/pq) = β , из которого можно найти

Nβ =

[

(xβ/2)2p(1-p)         Δ2

]

+ 1,

, где Φ₀(x_β/2) = β/2. Но величины p и q = 1 - p заранее неизвестны, поэтому заменим величину p(1 - p) на ее максимальное значение, которое, очевидно, достигается при p = 1/2 и равно 1/4. Поэтому выбираем N_β = [(x_β/2)²/(2Δ²)] + 1. В этом примере для вычисления N_β необходимо задать уровень доверительной вероятности β (вторичный по отношению к p) и точность Δ вычисления p. Например, при p > 0.8 обычно задают Δ = (1-p)/10 и β = 1-Δ.

Пример 2.   Отметим, что в предыдущем примере гарантирующее число испытаний было выбрано априори, т.е. до опыта, и одним и тем же для всех вычисляемых вероятностей p. Но интуитивно ясно, что N_β должно быть различным в зависимости от значения p и от реализации числа успешных испытаний M. Например, пусть в серии из N испытаний все испытания оказались успешными, т.е. M = N. Найдем для этого случая зависимость N_β(p) гарантирующего числа испытаний от доверительной вероятности β и вычисляемой вероятности p. С этой целью рассмотрим, как и в примере 1, нормированную частоту

* WN,

но потребуем выполнение другого вероятностного условия:

P{	* Wn ≤ xβ} = β,

где x_β - квантиль уровня β для стандартного нормального распределения N(0,1). Это условие отражает наше желание добиться того, чтобы частота W_n была не меньше оцениваемой вероятности p с некоторой точностью, т.е.

* WN

Δ = (WN - p)

√

N       p(1 - p)

= (

N-K   N

- p)

√

N       p(1 - p)

≤ xβ ,

где

K

Δ = N - M

есть число неуспешных испытаний. Действительно, это неравенство эквивалентно следующему:

WN ≤ p + xβ

√

p(1 - p)     N

.

Если последнее неравенство выполнено для конкретно реализовавшегося значения частоты, то можно надеяться (с доверительной вероятностью β), что полученное значение W_N лишь незначительно меньше вычисляемой вероятности p (см. более подробно о понятии доверительного интервала в Л15.Р4.О2). Выбирая минимальное N, удовлетворяющее полученному неравенству, т.е. решая соответствующее квадратное уравнение (N - K - Np)² = N(x_β)²p(1 - p), находим

Nβ(p,K) =

[

__________ 2K+(xβ)2p+xβ√4Kp+p2(xβ)2                 2(1 - p)

]

+ 1.

В частности, если все испытания оказались успешными (K = 0), то

Nβ(p,0) =

[

(xβ)2 p   1 - p

]

+ 1.

Например, если при вычислении вероятности p = 0.99 положить, как рекомендуется в примере 1, δ = 0.001, β = 0.999, то по формуле из этого примера получаем априорное значение гарантирующего числа испытаний N_β = 10 300 000, т.к. в этом случае x_β = 3.209. При тех же данных апостериорное значение гарантирующего числа непрерывно успешных испытаний оказывается намного меньше: N_β(p,0) = 1 030.

Лекция 13.
Оглавление