Введение в теорию вероятностей. Метод Монте-Карло. Сравнение оценок

Введение в теорию вероятностей Лекция 5 Метод Монте-Карло Метод Монте-Карло — метод решения математических задач при помощи случайных чисел. Знакомство с методом начнем с рассмотрения задачи численного интегрирования функции φ(x), заданной R1 на отрезке [0, 1]. Как вычислить приближенно интеграл I = 0 φ(x)dx? Пожалуй, простейший спо1 Pn соб — метод прямоугольников. Он состоит в замене I на интегральнуюvсумму In = n i=1 φ(xi ), где xi = i−1/2 — это «узлы» равномерной сетки, т. е. середины интервалов разбиения отрезка [0, 1] n на n равных частей. При условии, что φ(x) дважды непрерывно дифференцируема, можно показать, что погрешность метода прямоугольников δn = |I − In | оцениваться сверху так: δn ≤ M , 24n2 где M = max |φ00 (x)|. x∈[0,1] Таким образом, для гладких функций погрешность метода прямоугольников имеет порядок малости 1/n2 . Метод Монте-Карло для вычисления интеграла I отличается от метода прямоугольников тем, что в качестве «узлов» используются случайные числа y1 , . . . , yn . Пусть случайные величины X1 , . . . , Xn — независимы и равномерно распределены на [0, 1]. Тогда Z +∞ Eφ(Xi ) = −∞ Z 1 φ(u)du = I. φ(u)fXi (u)du = Согласно закону больших чисел n 1X P Iˆn = φ(Xi ) → I n при n → ∞. i=1 P Здесь → обозначает сходимость по вероятности: для любого ε > 0 P(|Iˆn − I| > ε) → 0 при n → ∞. Таким образом, с ростом n погрешность приближения должна стремиться к нулю. Кроме того, если σ 2 = Varφ(Xi ) < ∞, то по центральной предельной теореме (ЦПТ) ! √ ˆ n|In − I| P x≤ ≤ y → P(x ≤ Z ≤ y), σ где Z ∼ N (0, 1). Если положить x = −3 и y = 3, то мы получим P (−3 ≤ Z ≤ 3) ≈ 0.997. Мы сейчас получили, что при достаточно больших n (когда ЦПТ дает хорошую аппроксимацию) выполняется 1 Введение в теорию вероятностей Лекция 5 неравенство 3σ |Iˆn − I| ≤ √ n с вероятностью близкой к 1. √ Заметим, что ошибка оценки интеграла |Iˆn − I| убывает как 1/ n с ростом n. Говорят, что Iˆn является «оценкой» I, а интервал 3σ 3σ Iˆn − √ ≤ I ≤ Iˆn + √ n n называют «доверительным интервалом» оценки Iˆn . Отметим, что в выводе мы воспользовались так называемым «правилом трех сигм»: случайная величина X ∼ N (µ, σ 2 ) принимает значения из отрезка [µ − 3σ, µ + 3σ] с вероятностью 0.997, которую зачастую не отличают от 1. Неразумно использовать метод Монте-Карло для вычисления одномерных интегралов — для этого существуют квадратурные формулы, простейшая из которых — рассмотренная выше формула метода прямоугольников. Дело в том, что погрешность метода Монте-Карло больше, чем погрешность метода прямоугольников (да и верна эта оценка погрешности лишь с некоторой вероятностью). Тем не менее, метод Монте-Карло (или его модификации) часто оказывается единственным численным методом, позволяющим решить задачу вычисления интеграла большой кратности. Дело в том, что число «узлов» сетки возрастает как nk , где k — кратность интеграла (так называемое «проклятие размерности»). Так, чтобы найти интеграл по десятимерному кубу, используя в качестве «узлов» только его вершины, надо 210 = 1024 раза вычислить значение интегрируемой функции. В практических задачах эти вычисления могут оказаться довольно долгими, например, когда для расчета значений требуется численное решение систем нелинейных или дифференциальных уравнений. Напротив, метод Монте-Карло не зависит от размерности: чтобы найти приближенное значение интеграла Z 1 Z 1 ... φ(u1 , . . . , un )du1 . . . dun √ с точностью порядка 1/ n достаточно случайно набросать n точек в k-мерный единичный куб (разбив случайные числа на группы из k элементов) и вычислить среднее арифметическое значений φ в этих точках. В частности, если функция — индикатор некоторой области D, то с помощью метода Монте-Карло можно приближенно определить объем этой области. Например, частота случайных точек, попавших под дугу окружности будет служить приближением к π/4. Упражнение 1. Сколько случайных точек надо бросить в единичный квадрат, чтобы получить 2 Введение в теорию вероятностей Лекция 5 площадь под дугой окружности (см. рисунок выше) с точностью 0.001 и с вероятностью 0.997? Сравнение оценок Анализируемые методами математической статистики данные обычно рассматриваются как реализация выборки из некоторого распределения, известного с точностью до параметра (или нескольких параметров). При таком подходе для определения распределения, наиболее подходящего для описания данных, достаточно уметь оценивать значение параметра по реализации. В этой главе будет рассказано, как сравнивать различные оценки по точности. Статистическая модель Эксперимент. Пусть θ — некоторое неизвестное положительное число. Ниже приведены (с точностью до 0, 1) координаты xi десяти точек, взятых наудачу из отрезка [0, θ]. 3.5 3.2 25.6 8.8 11.6 26.6 18.2 0.4 12.3 30.1. Попробуйте угадать значение параметра θ, на котором изображены эти точки. С формальной точки зрения мы имеем дело со следующей моделью: набор xi — это реализация независимых и равномерно распределенных на отрезке [0, θ] случайных величин Xi c функцией распределения  0, если x ≤ 0,  Fθ (x) = x/θ, если 0 < x < θ,   1, если x ≥ θ. Здесь θ ∈ Θ = (0, +∞) — неизвестный параметр масштаба. Статистическая модель. В общем случае задается семейство функций распределения {Fθ (x), θ ∈ Θ}, где Θ — множество возможных значений параметра; данные x1 , . . . , xn рассматриваются как реализация выборки X1 , . . . , Xn , элементы которой имеют функцию распределения Fθ0 (x) при некотором неизвестном значении θ0 ∈ Θ. Задача состоит в том, чтобы оценить (восстановить) θ0 по выборке x1 , . . . , xn , по возможности, наиболее точно. Как «угадать» задуманное значение, основываясь на наблюдениях x1 , . . . , xn ? Будем оценивать θ0 при помощи некоторых функций θ̂ от n переменных x1 , . . . , xn , которые называются оценками или статистиками. Для приведенных выше данных эксперимента в качестве оценок неизвестного параметра масштаба можно использовать, скажем, θ̂1 = x(n) = max{x1 , . . . , xn } и θ̂2 = 2(x1 + . . . + xn )/n. Интуитивно понятно, что при увеличении n каждая из оценок будет приближаться именно к тому значению θ, с которым моделировалась выборка. Но какая из них точнее? Каким образом вообще можно сравнивать оценки? Прежде чем дать ответы на эти вопросы, познакомимся с важнейшими свойствами оценок — несмещенностью и состоятельностью. 3 Введение в теорию вероятностей Лекция 5 Несмещенность и состоятельность Определение. Оценка θ̂(x1 , . . . , xn ) параметра θ называется несмещенной, если Eθ θ̂(X1 , . . . , Xn ) = θ для всех θ ∈ Θ. Здесь индекс θ у Eθ означает, что имеется в виду математическое ожидание случайной величины θ̂(X1 , . . . , Xn ), где Xi распределены с функцией распределения Fθ (x). В дальнейшем этот индекс будет опускаться, чтобы формулы не выглядели слишком громоздко. Замечание. Важно, чтобы условие несмещенности выполнялось для всех θ ∈ Θ. Тривиальный контрпример: оценка θ̂(x1 , . . . , xn ) = 1, идеальная при θ = 1, при других значениях θ имеет смещение b(θ) = Eθ̂ − θ = 1 − θ. Иногда представляет интерес получение оценки не для самого параметра θ, а для некоторой заданной функции φ(θ). Пример 1. Для выборочного контроля из партии готовой продукции отобраны n приборов. Пусть величины X1 , . . . , Xn — их времена работы до поломки. Допустим, что Xi одинаково показательно распределены с неизвестным параметром θ : Fθ (x) = 1 − e−θx , x > 0. Требуется оценить среднее время до поломки прибора Z +∞ Z 1 +∞ −y 1 −θx φ(θ) = EX1 = θ xe dx = ye dy = . θ 0 θ По свойствам математического ожидания выборочное среднее X = оценкой для функции φ(θ): EX = φ(θ). X1 +...+Xn n будет несмещенной Пример 2. Рассмотрим выборку из какого-либо распределения с двумя параметрами µ и σ, где µ = EX1 и σ 2 = VarX1 (скажем, нормального закона N (µ, σ 2 )). По свойствам математического ожидания выборочное среднее X несмещенно оценивает параметр µ. В качестве оценки для неизвестной дисперсии φ(σ) = σ 2 можно взять выборочную дисперсию S2 = n n i=1 i=1 1X 1X 2 2 (Xi − X)2 = Xi − X . n n Однако, оценка S 2 имеет смещение. Действительно, так как случайные величины Xi независимы и одинаково распределены, то, применяя свойства математического ожидания, получаем:   n n n n n X X 1 1 1X 1 X 1 X 2 2 2 2   ES = E Xi − 2 Xi Xj = EXi − 2 EXi − 2 EXi EXj n n n n n i=1 = EX12 − i,j=1 i=1 i=1 i6=j 1 n−1 n−1 EX12 − (EX1 )2 = VarX1 . n n n Чтобы устранить смещение, достаточно домножить S 2 на n/(n − 1). Само по себе свойство несмещенности не достаточно для того, чтобы оценка хорошо приближала неизвестный параметр. Например, первый элемент X1 выборки из закона Бернулли служит несмещенной оценкой для θ: EX1 = 0 · (1 − θ) + 1 · θ = θ. Однако, его возможные значения 0 и 1 даже не принадлежат Θ = (0, 1). Необходимо, чтобы погрешность приближения стремилась к нулю с увеличением размера выборки. Это свойство в математической статистике называется состоятельностью. 4 Введение в теорию вероятностей Лекция 5 Определение. Оценка θ̂(x1 , . . . , xn ) параметра θ называется состоятельной, если для всех θ ∈ Θ последовательность P θ̂n = θ̂(X1 , . . . , Xn ) → θ при n → ∞. Состоятельность оценки (а точнее — последовательности оценок {θ̂n }) означает концентрацию вероятностной массы около истинного значения параметра с ростом размера выборки n. Как установить, будет ли данная оценка состоятельной? Обычно оказывается полезным один из следующих трех способов: 1. Иногда удается доказать состоятельность, непосредственно вычисляя функцию распределения оценки. 2. Другой способ проверки состоит в использовании закона больших чисел и свойства сходимости: P если случайные величины ξn сходятся по вероятности к случайной величине ξ, то есть ξn → ξ, P и функция φ(x) непрерывна, то φ(ξn ) → φ(ξ). 3. Часто установить состоятельность помогает следующая лемма. Лемма. Если оценка θ̂n не смещена, Eθ̂n = θ, и дисперсия Varθ̂n стремятся к нулю при n → ∞, то оценка θ̂n состоятельна. Доказательство. Согласно неравенству Чебышева: P(|θ̂n − θ| > ε) ≤ Varθ̂n → 0 при n → ∞. ε2 Лемма доказана. Упражнение 2. Для случайных величин Xi , i = 1, . . . , n, взятых наудачу из отрезка [0, θ], докажите состоятельность оценки θ̂ = X(n) = max{X1 , . . . , Xn } первым способом. Упражнение 3. Для случайных величин, распределенных экспоненциально с неизвестным параметром θ > 0 (см. Пример 1), докажите состоятельность оценки θ̂ = 1/X параметра θ вторым способом. Упражнение 4. Для случайных величин Xi , i = 1, . . . , n, взятых из распределения N (θ, 1), докажите состоятельность оценки θ̂ = X = (X1 + . . . + Xn )/n третьим способом. Методы получения оценок В этой главе рассматриваются несколько методов получения оценок параметров статистических моделей, в том числе — метод моментов и метод максимального правдоподобия. Plug-in оценки Если нам необходимо оценить параметр θ и θ = Eφ(X1 ) для некоторой функции φ(x), то мы можем рассмотреть оценку вида n 1X θ̂ = φ(Xi ). n i=1 Нам уже встречались оценки такого типа. Например, в модели Бернулли, когда Xi имею вероятность успеха с неизвестным параметром θ ∈ (0, 1), мы оценивали θ = EX1 с помощью среднего X. По закону больших чисел такие оценки оказываются несмещенными и состоятельными. Если Varφ(X1 ) < +∞, то можно построить доверительный интервал для plug-in оценки (как мы делали в главе про Монте-Карло). 5 Введение в теорию вероятностей Лекция 5 Метод моментов Моментом k-го порядка случайной величины X называется величина αk = EX k . Моменты существуют не всегда. Например, у закона Коши математическое ожидание α1 не определено. n 1X k P Xi . Если момент αk существует, то в силу закона больших чисел Ak → αk (это n i=1 plug-in оценки моментов). Поэтому для реализации x1 , . . . , xn выборки достаточно большого размеn 1X k ра можно утверждать, что ak = xi ≈ αk , т. е. эмпирические моменты k-го порядка ak близки n i=1 к теоретическим моментам αk . На этом соображении основывается так называемый метод моментов. Положим Ak = Допустим, что распределение элементов выборки зависит от m неизвестных параметров θ1 , . . . , θm , где вектор θ = (θ1 , . . . , θm ) принадлежит некоторой области Θ в Rm . Пусть E|X|m < ∞ для всех θ ∈ Θ (отсюда следует конечность всех моментов до m из неравенства Ляпунова). Тогда существуют все αk = αk (θ), k = 1, . . . , m, и можно записать систему из m (вообще говоря, нелинейных) уравнений αk (θ) = ak , k = 1, . . . , m. Предположим, что левая часть системы задает взаимно однозначное отображение g : Θ → B, где B — некоторая область в Rm , и что обратное отображение g −1 : B → Θ непрерывно. Другими словами, для всех (y1 , ..., ym ) из B система имеет единственное решение, которое непрерывно зависит от правой части. Компоненты решения θ̂ = (θ̂1 , ..., θ̂m ) при yk = Ak называются оценками метода моментов. Пример 3. Рассмотрим модель сдвига показательного закона, в которой плотностью распределения величин Xi служит функция pθ (x) = e−(x−θ) · 1{x≥θ} . Здесь Z α1 (θ) = EX1 = +∞ xe −(x−θ) Z dx = θ +∞ (y + θ)e−y dy = 1 + θ. Из уравнения 1 + θ = A1 = X, находим по методу моментов оценку θ̂ = X − 1. Какими статистическими свойствами обладают оценки, полученные методом моментов? Их состоятельность вытекает из непрерывности определенного выше отображения g −1 . Метод максимального правдоподобия Метод получил распространение после появления в 1912 г. статьи Р. Фишера, где было доказано, что получаемые этим методом оценки являются асимптотически наиболее точными при выполнении некоторых условий регулярности модели. Для знакомства с методом предположим для простоты, что элементы выборки Xi имеют дискретное распределение: f (x, θ) = P (X1 = x) (здесь θ = (θ1 , . . . , θm ) — вектор неизвестных параметров 6 Введение в теорию вероятностей Лекция 5 модели). Тогда совместная вероятность выборки f (x, θ) = f (x1 , θ) · . . . · f (xn , θ) зависит от n+m аргументов (здесь x = (x1 , ..., xn )). Рассматриваемая как функция от θ1 , . . . , θm при фиксированных значениях элементов выборки x1 , . . . , xn , она называется функцией правдоподобия и обычно обозначается через L(θ). Величину L(θ) можно считать мерой правдоподобия значения L(θ) при заданной реализации x. Представляется разумным в качестве оценок параметров θ1 , . . . , θm взять наиболее правдоподобные значения θ̃1 , . . . , θ̃m которые получаются при максимизации функции L(θ). Такие оценки называются оценками максимального правдоподобия (ОМП). Часто проще искать точку максимума функции ln L(θ), которая совпадает с θ̃ = (θ̃1 , . . . , θ̃m ) в силу монотонности логарифма. Пример 4. Для схемы Бернулли X1 , . . . , Xn с вероятностью «успеха» θ имеем: f (x, θ) = P (X1 = x) = θx (1 − θ)1−x , где x принимает значения 0 или 1. Поэтому функция правдоподобия L(θ) = θsn (1 − θ)n−sn , где sn = x1 + . . . + xn , представляет собой многочлен n-й степени. Найдем точку максимума ln L(θ) = sn ln θ + (n − sn ) ln(1 − θ). Дифференцируя по θ, получаем уравнение sn /θ − (n − sn )/(1 − θ) = 0, откуда θ̃ = sn /n = x̄. Таким образом, ОМП в схеме Бернулли — это частота «успехов» в реализации x1 , . . . , xn . В случае непрерывных моделей будем использовать обозначение f (x, θ) для плотности распределения случайной величины X1 . Пример 5. Рассмотрим модель сдвига показательного закона с плотностью f (x, θ) = e−(x−θ) 1{x≥θ} . В этом случае функция правдоподобия равна L(θ) = n Y f (xi , θ) = e−(x1 +...+xn ) enθ 1{x(1) ≥θ} . i=1 Отсюда (см. рис. ниже) получаем в качестве ОМП θ̃ = x(1) , которая отлична от оценки метода моментов θ̂ = x̄ − 1, найденной ранее для этой модели в Примере 3. Заметим также, что здесь L(θ) не является гладкой функцией, и поэтому ОМП нельзя вычислять, приравнивая нулю производную функции правдоподобия. В случае, когда L(θ) гладко зависит от θ1 , . . . , θm , оценки максимального правдоподобия являются 7 Введение в теорию вероятностей компонентами решения (вообще говоря, нелинейной) системы уравнений: n X ∂ ∂ ln L(θ) = ln f (xi , θ) = 0, ∂θj ∂θj j = 1, . . . , m. i=1 Список литературы [1] М. Б. Лагутин. Наглядная математическая статистика. Бином, 2009. 8 Лекция 5