W-функция Ламберта, ее приближения, и как она помогает строить доверительные интервалы

Читателям, которые не знают, что такое W-функция Ламберта, какие у нее есть ветви и как с ее помощью можно решать различные полиномиально-экспоненциальные уравнения, предлагаю ознакомиться с первой статьей.

В этой же статье мы сначала поговорим о том, как можно приближать значения различных ветвей W-функции Ламберта, а после применим полученные результаты для решения задачи поиска кратчайших доверительных интервалов. Немного о пороге входа: читателю понадобятся базовые знания в теории вероятностей, математической статистике и, что менее банально, комплексном анализе (хотя в более-менее серьезном виде он понадобится для получения только одного выражения, которое можно принять на веру).

Начнем с первого, что приходит на ум, когда говорят о приближениях - рядах Тейлора. Изучим вопрос того, как хорошо или плохо раскладывается W-функция в степенной ряд в окрестности нуля. Для этого сразу предлагается сбежать в комплексный анализ и ввести следующую функцию:

Теперь воспользуемся теоремой Лагранжа об обращении рядов. Приведем ее формулировку. Здесь не стоит пугаться того, что она дана в терминах ТФКП, просто с обращением функций и разложением их в ряд в этой науке проще, чем в действительном братце.

Теорема (Лагранжа, об обращении рядов). Пусть функция голоморфна в точке , причем , тогда в некоторой окрестности точки функция, обратная к может быть представлена в виде следующего ряда:

Посмотрим на то, как перепишется утверждение данной теоремы применительно к W-функции Ламберта:

Осталось найти соответствующие коэффициенты:

Как итог, получаем следующую симпатичную формулу, правда, с пока неизвестной областью применимости:

Найдем радиус сходимости этого степенного ряда с помощью теоремы Коши-Адамара и формулы Стирлинга:

Итак, полученный степенной ряд сходится в открытом круге .

Далее обсудим то, как можно использовать схожесть графика функции в окрестности точки с параболой. Для этого воспользуемся теоремой Тейлора с остаточным членом в форме Пеано с центром в точке -1:

Теперь так как нас интересуют обратные функции отражаем картинку относительно и получаем следующие предельные равенства:

Полученные приближения хороши при "не сильно больших" . Казалось бы, зачем они нам при наличии сходящегося ряда, однако ближе к граничным точкам интервала этот ряд сходится медленно и даже если взять 10 первых слагаемых, то около приближение будет хуже, чем при описанном выше методе, который к тому же оказывается полезным при оценке значений "минус первой" ветви.

Замечание. Можно раскладывать в Тейлора до члена третьей степени и дальше решать кубическое уравнение и получать лучшую аппроксимацию, но мы не будем этого делать, чтобы не устраивать радикальное безумие имени Кардано.

Теперь же мы приступаем к самой интересной с математической точки зрения части статьи. Если Вам не интересен комплексный анализ или не хочется вникать в математические выкладки (хотя дальше таковых будет не очень много), то сразу переходите на абзац с полученным результатом.

Итак, введем следующую функцию: . Она целая (то есть голоморфна во всей комплексной плоскости) и зависит от двух комплексных параметров и .

Утверждение 1. Существует такое число , что при функция имеет в точности один нуль внутри открытого круга .

Доказательство. Пусть на окружности . Заметим, что в силу того, что - это строго положительная непрерывная на компакте функция. Внутри этого круга ее единственный ноль -. Теперь явно укажем значение и покажем, что оно действительно подходит:

Далее нам понадобится теорема Руше, сформулируем ее.

Теорема (Руше). Пусть - ограниченная область с кусочно-гладкой границей , функции голоморфны в области , компактно содержащей и всюду на верно, что , тогда в области количество нулей с учетом кратности функции совпадает с количеством нулей с учетом кратности функции .

Звучит страшновато, но смысл простой: если на границе области одна из функций по модулю строго меньше другой, то и внутри области ее прибавление не оказывает достаточного влияния на другую функцию в смысле изменения количества нулей.

Вернемся к доказательству, которое по указанной выше теореме завершается следующей цепочкой неравенств:

Утверждение 2. Пусть - единственный корень в круге . Тогда верно следующее равенство:

Доказательство. Учитывая единственность особенности подыинтегральной функции, по теореме Коши о вычетах достаточно показать, что

Для этого рассмотрим предел, который равен "минус первому" члену в ряде Лорана. Преобразуем его, используя, что .

На этом доказательство утверждения завершено.

Сделаем следующий шаг: разложим знаменатель в ряд Тейлора, который в силу утверждения 1 будет абсолютно сходящимся (это позволяет не просто утверждать, что такое представление справедливо, но еще и менять порядок суммирования).

Меняем порядок суммирования с "зеленых" направлений на "синие" (см. рисунок)

Полученное в виде двойного ряда представление подставляем в изначальное интегральное равенство из утверждения 2 и меняем порядок операций.

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

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

Пока что это были отвлеченные от W-функции Ламберта рассуждения, но зато теперь мы готовы их с ней связать. Сделаем в основном Ламбертовом тождестве замену (берем основную ветвь комплексного логарифма): . Наводящие на именно такое преобразование соображения приведены в [1] вместе с англоязычным доказательством.

Сокращаем на и приходим к следующему равенству:

Очень похоже на ту функцию, с которой мы начали разложение, правда? Если теперь сделать замену, то мы получим ее в явном виде:

Для доказательства прекрасной формулы осталось заметить, что при достаточно больших и, наоборот, достаточно малых по модулю соответствующие значения и "зайдут" в границы, описанные в утверждении 1, а значит будет справедливо представление в виде двойного ряда, которое в терминах W-функции Ламберта перепишется следующим образом:

Кроме того, аналогичное представление верно и для "минус первой" ветви, только взять нужно и соответственно.

Замечание. На самом деле таким методом можно получить представление для любой ветви W-функции Ламберта (просто нас интересуют только "нулевая" и "минус первая"), нужно лишь брать другую ветвь внутреннего логарифма.

Снова вернемся в действительный мир и распишем несколько первых членов (в смысле степени знаменателя) ряда:

Подведем промежуточные итоги. Что мы умеем?

• приближать "нулевую" ветвь около 0

• приближать "нулевую" и "минус первую" ветвь около

• приближать "минус первую" ветвь всюду

Теперь перейдем к применению полученных результатов.

Раз мы заглянули в комплексный анализ, то нечего нам уже останавливаться и будем пользоваться удобным аппаратом характеристических функций.

Определение. Характеристической функцией случайной величины называется функция , определяемая следующий выражением:

Основная прелесть характеристических функций заключается в следующей теореме:

Теорема (Леви, о непрерывности). Пусть - последовательность случайных величин, тогда:

1. если по распределению при , то поточечно

2. если поточечно, причем непрерывна в нуле, то существует случайная величина такая, что и по распределению при

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

Определение. Говорят, что случайная величина имеет гамма-распределение с параметром масштаба и параметром формы , если ее плотность вероятности имеет следующий вид:

Обозначение:

Замечание. Как в википедии, так и в нескольких математических пакетах используется определение, получающееся из приведенного выше заменой, однако автор считает идейно правильным умножать на параметр масштаба, а не делить, хотя это вопроса вкуса и простоты выкладок в конкретной задаче.

Гамма-распределение прекрасно тем, что многие другие распределения являются его частными случаями, например, экспоненциальное распределение с параметром это .

Характеристическая функция случайной величины имеет вид:

Вывод этого факта оставляем читателю в качестве упражнения. Теперь мы можем сказать, например, как распределена случайная величина :

Утверждение 3. Пусть - независимые одинаково распределенные случайные величины, причем , тогда случайная величина .

Доказательство. Воспользуемся одним из свойств характеристических функций. Оно звучит так: характеристическая функция суммы независимых случайных величин равна произведению их характеристических функций. Это почти очевидное следствие свойств математических ожиданий и показательной функции.

Осталось заметить, что о распределении умноженной на число случайной величины мы уже говорили.

Для начала освежим в памяти понятие доверительного интервала. Пусть выборка состоит из независимых одинаково распределенных случайных величин , а множество распределений параметризовано с помощью .

Определение. Интервал называется доверительным интервалом с уровнем доверия , если для любого вероятность того, что он накроет неизвестный параметр, равна :

Постановка задачи. Пусть и дополнительно . Нужно построить кратчайший доверительный интервал для параметра .

Решение. Сначала разберемся с тем, что доверительный интервал должен быть кратчайшим, то есть нужно решить следующую оптимизационную задачу (здесь - плотность распределения случайной величины, как-то слепленной из выборки, пока это не важно, нужно уточнить другой момент):

Будем решать ее методом множителей Лагранжа, составим функционал и будем его безусловно минимизировать:

Необходимое условие экстремума запишется следующим образом:

То есть если минимум достигается, то значения плотности в концах доверительного интервала должны быть равны. Запомним это.

Решать поставленную задачу будем с помощью метода центральной случайной величины: для этого нужно найти такую случайную величину , распределение которой не зависит от параметра , тогда можно, используя, например, квантили этого распределения и монотонную зависимость от параметра, перейти от интервала для центральной случайной величины к интервалу для параметра. В нашей задаче в качестве центральной случайной величины будем использовать (вид распределения следует из утверждения 3):

Тогда верны следующие равенства:

Вспоминаем рассуждения в задаче минимизации и получаем переформулированную постановку: нужно найти такой уровень , что если - упорядоченные решения уравнения , то

Делать это можно, например, с помощью бинарного поиска с заданной точностью, начиная из и .

Замечание. Дополнительное условие на величину было поставлено, чтобы плотность имела именно такой вид, как на рисунке, или, строго говоря, не имела особенностей в нуле.

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

Обратим внимание на то, что , поэтому при реализации алгоритма мы можем пользоваться приближениями, полученными в первой части статьи.

Теперь у нас есть два варианта, как считать определенный интеграл:

1. Пользоваться численными методами (прямоугольников, трапеций, парабол) и получать погрешность, во-первых, из-за границ интегрирования, во-вторых, из-за самих численных методов.

2. Предположить, что , явно выписать первообразную и вычислить нужную величину по формуле Ньютона-Лейбница. Тогда погрешность будет накапливаться из-за возведения значений границ интегрирования в степени.

С первым методом все понятно, рассмотрим подробнее второй.

Утверждение 4.

Доказательство. Оставляем читателю в качестве упражнения по методу математической индукции и интегрированию по частям.

Вычисляем определенный интеграл:

Далее пользуемся основным Ламбертовым тождеством:

Делая замену , получаем чуть более простой вид:

Заметим, что множитель перед суммой симпатично выглядит, если раскрыть его по определению:

Подставляем и делаем замену индекса суммирования:

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

Теорема. Пусть - четное число, - интервал для , найденный описанным выше методом, тогда для любого

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

В [2] можно найти асимптотическую формулу для квантилей, которая справедлива в условиях теоремы. Применяем ее и получаем, что разность в правой части неравенства есть , что в свою очередь позволяет расписать эту вероятность по определению и оценить сверху сходящейся к 0 последовательностью.

Итак, перейдем к результатам тестирования. Метод с явным вычислением первообразной хорош в тех случаях, когда - невелико, в противном случае его начинает по понятным причинам штормить вплоть до удивительных явлений в виде вероятности, большей единицы. Далее приводятся результаты экспериментов для метода, в котором при подсчете интеграла использовались численные методы, а именно функция . Всюду в экспериментах проводилось 2000 запусков, уровень доверия .

Раз-р выборки,

Ср. длина инт-ла

0.5

0.1

50

1.03

100

0.67

200

0.47

1

0.2

50

1.34

100

0.92

200

0.65

2

0.5

50

1.64

100

1.15

200

0.8

3

0.75

50

2.0

100

1.39

150

1.13

0.25

1

25

0.21

50

0.14

100

0.1

1. Corless, R.M., Gonnet, G.H., Hare, D.E.G. et al. On the LambertW function. Adv Comput Math 5, 329–359 (1996)

2. В. И. Пагурова, “О вычислении квантилей Γ-распределения”, Теория вероятн. и ее примен., 10:4 (1965), 746–749; Theory Probab. Appl., 10:4 (1965), 677–680