Хабрахабр

В трёх статьях о наименьших квадратах: ликбез по теории вероятностей

Полтора года назад я опубликовал статью «Математика на пальцах: методы наименьших квадратов», которая получила весьма приличный отклик, который, в том числе, заключался в том, что я предложил нарисовать сову. Ну, раз сова, значит, нужно объяснять ещё раз. Через неделю ровно на эту тему я начну читать несколько лекций студентам-геологам; пользуюсь случаем, излагаю тут (адаптированные) основные тезисы в качестве черновика. Моей основной целью не является дать готовый рецепт из книги о вкусной и здоровой пищи, но рассказать, почему он таков и что ещё находится в соответствующем разделе, ведь связи между разными разделами математики — это самое интересное!

На данный момент я предполагаю разбить текст на три статьи:

  • 1. Ликбез по теории вероятностей и как она связана с методами наименьших квадратов
  • 2. Наименьшие квадраты, простейший случай, и как их программировать
  • 3. Нелинейные задачи

Я зайду к наименьшим квадратам чуть сбоку, через принцип максимума правдоподобия, а он требует минимального ориентирования в теории вероятностей. Данный текст рассчитан на третий курс нашего факультета геологии, что означает, (с точки зрения задействованного матаппарата!) что заинтересованный старшеклассник при соответствующем усердии должен суметь в нём разобраться.
Однажды мне задали вопрос, верю ли я в теорию эволюции. Прямо сейчас сделайте паузу, подумайте, как вы на него ответите.

Научная теория имеет мало общего с верой.
Лично я, опешив, ответил, что нахожу её правдоподобной, и что вопрос веры тут вообще не встаёт. Более того, критерий Поппера требует от научной теории иметь возможность опровержения. Если кратко, то теория лишь строит модель окружающего нас мира, нет необходимости в неё верить. Например, если вы генетически модифицируете сельскохозяйственные культуры таким образом, что они сами будут производить пестициды, то вполне логично, что будут появляться устойчивые к ним насекомые. А ещё состоятельная теория должна обладать, в первую очередь, предсказательной силой. Основываясь на теории эволюции, соответствующее моделирование сделало такой прогноз, и он, похоже, подтвердился. Однако существенно менее очевидно, что этот процесс может быть замедлен благодаря выращиванию обычных растений бок о бок с генномодифицированными.

А при чём тут наименьшие квадраты?

Как я упомянул ранее, я зайду к наименьшим квадратам через принцип максимального правдоподобия. Давайте проиллюстрируем на примере. Предположим, что нам интересны данные о росте пингвинов, но мы можем измерить лишь несколько этих прекрасных птиц. Вполне логично ввести в задачу модель распределения роста — чаще всего она нормальная. Нормальное распределение характеризуется двумя параметрами — средним значением и среднеквадратичным отклонением. Для каждого фиксированного значения параметров мы можем посчитать вероятность того, что будут сгенерированы именно те измерения, что мы произвели. Дальше, варьируя параметры, найдём те, что максимизируют вероятность.

Чуть ниже мы «на пальцах» определим понятие вероятности и правдоподобия, но сначала я хотел бы заострить внимание на другом аспекте. Таким образом, для работы с максимальным правдоподобием нам нужно оперировать в понятиях теории вероятностей. Я на удивление редко вижу людей, которые задумываются над словом «теория» в словосочетании «теория вероятностей».

Что изучает теорвер?

По поводу того, каковы истоки, значения и область применимости вероятностных оценок, уже не первую сотню лет идут яростные споры. Например, Брюно Де Финетти заявил, что вероятность есть не что иное, как субъективный анализ вероятности того, что что-то произойдёт, и что эта вероятность не существует вне ума. Это готовность человека делать ставки на что-то происходящее. Это мнение прямо противоположно взгляду классиков / фреквентистов на вероятность конкретного результата события, в котором предполагается, что одно и то же событие может повторяться многократно, а «вероятность» конкретного результата связана с частотой выпадения конкретного результата во время повторных испытаний. Помимо субъективистов и фреквентистов есть ещё и объективисты, утверждающие, что вероятности — это реальные аспекты универсума, а не просто описания степени уверенности наблюдателя.

Давайте приведём косвенный аргумент, с субъективистской точки зрения, в пользу теории вероятностей, построенной на аксиомах Колмогорова. Как бы то ни было, но все три научные школы на практике используют один и тот же аппарат, основанный на аксиомах Колмогорова. Пусть у нас будут два события: a = чемпионом станет команда Уругвая, b = чемпионом станет команда Германии. Приведём сами аксиомы чуть позже, а для начала предположим, что у нас есть букмекер, который принимает ставки на следующий чемпионат мира по футболу. Очевидно, что и Германия, и Уругвай победить одновременно не могут, поэтому шанс a∧b нулевой. Букмекер оценивает шансы команды Уругвая победить в 40%, шансы команды Германии в 30%. Давайте это запишем в следующем виде: Ну и заодно букмекер считает, что вероятность того, что победит либо Уругвай, либо Германия (а не Аргентина или Австралия) составляет 80%.

4, то есть, P(a) = 0. Если букмекер утверждает, что его степень уверенности в событии а равняется 0. Это значит, что игрок может сделать ставку на то, что событие а произойдёт, поставив четыре рубля против шести рублей букмекера. 4, то игрок может выбрать, будет ли он делать ставку за или против высказывания а, ставя суммы, которые совместимы со степенью уверенности букмекера. Или же игрок может поставить шесть рублей вместо четырёх рублей букмекера на то, что событие а не произойдёт.

Более того, в данном конкретном примере у игрока есть стратегия, при которой букмекер всегда теряет деньги. Если степень уверенности букмекера неточно отражает состояние мира, то можно рассчитывать на то, что в долговременной перспективе он будет проигрывать деньги игрокам, чьи убеждения более точны. Давайте её проиллюстрируем:

Обратите внимание, что в рассмотрение суммы выигрыша в принципе не входит то, являются ли Уругвай или Германия фаворитами чемпионата, проигрыш букмекера обеспечен! Игрок делает три ставки, и какой бы ни был исход чемпионата, он всегда в выигрыше. К этой ситуации привело то, что букмекер не руководствовался азами теории вероятностей, нарушив третью аксиому Колмогорова, давайте приведём их все три:

В текстовом виде они выглядят так:

  • 1. Все вероятности находятся в пределах от 0 до 1
  • 2. Безусловно истинные высказывания имеют вероятность 1, а безусловно ложные вероятность 0.
  • 3. Третья аксиома — это аксиома дизъюнкции, её легко интуитивно понять, отметив, что те случаи, когда высказывание a является истинным, вместе с теми случаями, когда b является истинным, безусловно охватывает все те случаи, когда истинно высказывание a∨b; но в сумме двух множеств случаев их пересечение встречается дважды, поэтому небходимо вычесть P(a∧b).

В 1931 году де Финетти доказал очень сильное утверждение:

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

Аксиомы вероятностей могут рассматриваться как ограничивающее множество вероятностных убеждений, которых может придерживаться некоторый агент. Обратите внимание, что из следования букмекера аксиомам Колмогорова не вытекает то, что он будет выигрывать (оставим в стороне вопросы комиссионных), но если им не следовать, то он будет гарантированно проигрывать. Отметим, что в пользу применения вероятностей были выдвинуты и другие аргументы; но именно практический успех систем формирования рассуждений, основанных на теории вероятностей, оказался привлекательным стимулом, вызвавшим пересмотр многих взглядов.

Весь теорвер построен лишь на трёх аксиомах; во всех трёх участвует некая магическая функция P. Итак, мы чуть-чуть приоткрыли завесу того, почему теорвер может иметь смысл, но какими именно объектами он манипулирует? Давайте попробуем посмотреть, работает ли площадь для определения вероятности. Более того, глядя на эти аксиомы, мне это очень напоминает функцию площади фигуры.

Определим слово «вероятность события» как «площадь соответствующего подмножества». Определим слово «событие» как «подмножество единичного квадрата». Шансы того, что пуля попадёт в данное множество, прямо пропорциональны площади множества. Грубо говоря, у нас есть большая картонная мишень, и мы, закрыв глаза, в неё стреляем. Из нашего определения вероятности следует, что идеально попасть в точку невозможно (наша пуля — это материальная точка). Достоверное событие в данном случае — весь квадрат, а заведомо ложное, например, любая точка квадрата. Давайте проиллюстрируем все три аксиомы: Я очень люблю картинки, и рисую их много, и теорвер не является исключением!

Достоверное событие — это весь квадрат, а заведомо ложное — любое множество нулевой площади. Итак, первая аксиома выполнена: площадь неотрицательна, и не может превысить единицы. И с дизъюнкицей работает отлично!

Пример первый: подбрасывание монетки

Давайте рассмотрим простейший пример с подбрасыванием монеты, он же схема Бернулли. Проводится n опытов, в каждом из которых может произойти одно из двух событий («успех» или «неудача»), одно с вероятностью p, второе с вероятностью 1-p. Наша задача состоит в том, чтобы найти вероятность получения ровно k успехов в этих n опытах. Эту вероятность нам даёт формула Бернулли:

5), подбрасываем её десять раз (n=10), и считаем, сколько раз выпадает решка: Возьмём обычную монету (p=0.

Вот так выглядит график плотности вероятности:

5), и также зафиксировали количество экспериментов (10), то возможное количество «успехов» может быть любым целым числом между 0 и 10, однако эти исходы не являются равновероятными. Таким образом, если мы зафиксировали вероятность наступления «успеха» (0. Например, вероятность насчитать семь решек составляет примерно 12%. Вполне очевидно, что получить пять «успехов» гораздо вероятнее, нежели ни одного.

У нас есть реальная монетка, но её распределение априорной вероятности «успех»/«неудача» мы не знаем. А теперь давайте посмотрим на эту же задачу с другой стороны. Например, у нас выпало семь решек. Однако мы можем её подкинуть десять раз и посчитать количество «успехов». Как нам это поможет оценить p?

Мы можем попробовать зафиксировать в формуле Бернулли n=10 и k=7, оставив p свободным параметром:

Я даже буковку сменил у функции, теперь это L (от англ. Тогда формула Бернулли может быть интерпретирована как правдоподобие оцениваемого параметра, (в данном случае p). То есть, правдоподобие — это вероятность сгенерировать данные наблюдения (7 решек из 10 опытов) для заданного значения параметра(-ов). likehood).

5) при условии выпадения семи решек из десяти бросков примерно равняется 12%. Например, правдоподобие сбалансированной монеты (p=0. Можно построить график функции L:

В данном конкретном случае у нас функция одной переменной, мы ищем её максимум. Итак, мы ищем такое значение параметров, которое максимизирует правдоподобие получения тех наблюдений, что у нас есть. Логарифм — функция строго монотонная, так что максимизация одного и другого — это строго одно и то же. Для того, чтобы было проще искать, я буду искать максимум не L, а log L. Итак, мы ищем максимум вот этой функции: А логарифм нам разбивает произведение на сумму, которую сильно удобнее дифференцировать.

Для этого приравняем нулю её производную:

Производная log x = 1/x, получаем:

То есть, максимум правдоподобия (примерно 27%) достигается в точке

На всякий случай посчитаем вторую производную:

7 она отрицательна, так что эта точка действительно является максимумом функции L. В точке p=0.

7: А вот так выглядит плотность вероятности для схемы Бернулли с p = 0.

Пример второй: АЦП

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

Какое из них брать? Например, если мы измерим напряжение на батарейке N раз, то получим N разных измерений. Итак, пусть у нас будет N величин Uj: Все!

Вот так выглядит плотность вероятности: Предположим, что каждое измерение Uj равно идеальной величине, плюс гауссов шум, который характеризуется двумя параметрами — положение гауссова колокола и его «ширина».

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

Ну а дальше всё строго как раньше, приравняем нулю частные производные по параметрам, которые мы ищем:

Получим, что наиболее вероятная оценка неизвестной величины U может быть найдена как среднее от всех измерений:

Ну а наиболее вероятный параметр сигма это обычное среднеквадратичное отклонение:

На мой вкус, стоило. Стоило ли заморачиваться, чтобы получить в ответе простое среднее от всех измерений? Например, усреднение АЦП. Кстати, усреднение многократных измерений постоянной величины с целью увеличения точности измерений — это вполне стандартная практика. К слову сказать, для этого гауссовости шума не нужно, достаточно того, что шум несмешённый.

Пример третий, и снова одномерный

Продолжаем разговор, давайте возьмём тот же самый пример, но чуть его усложним. Мы хотим измерить сопротивление некоего резистора. При помощи лабораторного блока питания мы умеем пропустить через него некоторое эталонное количество ампер, и измерить напряжение, которое для этого понадобится. То есть, у нас на вход нашего оценщика сопротивления будет N пар чисел (Ij, Uj).

Нарисуем эти точки на графике; закон Ома нам говорит о том, что мы ищем наклон синей прямой.

Запишем выражение для правдоподобия параметра R:

И опять приравняем нулю соответствующую частную производную:

Тогда наиболее правдоподобное сопротивление R может быть найдено по следующей формуле:

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

Наверняка вы уже обратили внимание, что в последних двух примерах максимизация логарифма правдоподобия эквивалентна минимизации суммы квадратов ошибки оценивания. Давайте рассмотрим ещё один пример. Возьмём калибровку безмена при помощи эталонных грузов. Пусть у нас есть N эталонных грузов массой xj, повесим их на безмен и измерим длину пружины, получим N длин пружины yj:

Пусть жёсткость пружины — это параметр a, ну а растяжение пружины под собственным весом — это параметр b. Закон Гука нам говорит, что растяжение пружины линейно зависит от приложенной силы, а в эту силу у нас входит вес грузов и вес самой пружины. Тогда мы можем вот таким образом записать выражение правдоподобия наших измерений (по-прежнему под гипотезой гауссового шума измерений):

Максимизация правдоподобия L эквивалентна минимизации суммы квадратов ошибок оценивания, то есть, мы можем искать минимум функции S, определённой следующим образом:

Иными словами, мы ищем такую прямую, которая минимизирует сумму квадратов длин зелёных отрезков:

Ну а дальше никаких сюрпризов, приравниваем к нулю частные производные:

Получаем систему из двух линейных уравнений с двумя неизвестными:

Вспоминаем седьмой класс школы и выписываем решение:

Методы наименьших квадратов являются частным случаем максимизации правдоподобия для тех случаев, когда плотность вероятности является гауссовой. В случае же когда плотность (совсем) не гауссова, МНК дают оценку, отличающуюся от оценки MLE (maximum likehood estimation). Кстати, в своё время Гаусс выдвигал гипотезу, что распределение не играет роли, важна только независимость испытаний.

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

Теги
Показать больше

Похожие статьи

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»
Закрыть