Карта → событие

Париж 1805; Гёттинген, 1809

Лежандр и Гаусс: метод наименьших квадратов

Теория вероятностей

Задача

Астроном делает $n$ измерений, из которых надо определить $m$ параметров орбиты, причём $n>m$. Система уравнений переопределена: точного решения нет, потому что измерения содержат ошибки.

До 1805 года действовали кустарно: отбрасывали «плохие» наблюдения, группировали уравнения и складывали, усредняли по подгруппам. Каждый астроном имел свой рецепт, обоснования не было ни у кого.

Решение Лежандра

Адриен Мари Лежандр. Единственное достоверное изображение — карикатура из альбома Института Франции
Адриен Мари Лежандр. Единственное достоверное изображение — карикатура из альбома Института Францииавтор неизвестен · Public domain

Адриен Мари Лежандр, приложение к «Nouvelles méthodes pour la détermination des orbites des comètes» (Париж, 1805). Девять страниц под заголовком «Sur la méthode des moindres quarrés» — он же дал методу имя.

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

$$\Phi(\theta) = \sum_{i=1}^{n} r_{i}^{2} = \sum_{i=1}^{n}\left(y_{i} - f(x_{i};\theta)\right)^{2} \longrightarrow \min.$$

Почему квадраты, а не модули. Лежандр приводит практические доводы: квадраты дифференцируемы (модуль в нуле — нет), поэтому условие минимума даёт линейную систему, решаемую до конца. Кроме того, квадраты сильнее штрафуют крупные промахи.

Нормальные уравнения. Для линейной модели $y = X\theta + \varepsilon$ условие $\nabla\Phi = 0$ даёт

$$X^{\top}X\,\hat\theta = X^{\top}y \quad\Longrightarrow\quad \hat\theta = \left(X^{\top}X\right)^{-1}X^{\top}y.$$

Простейший случай, который стоит проделать на уроке. Модель $y=\theta$ (просто оцениваем неизвестную величину): $\Phi(\theta)=\sum(y_i-\theta)^2$, производная $-2\sum(y_i-\theta)=0$, откуда

$$\hat\theta = \frac{1}{n}\sum y_{i} = \bar y.$$

Метод наименьших квадратов, применённый к одной величине, даёт среднее арифметическое. Запомним это — сейчас пригодится.

Приоритетный спор

Карл Фридрих Гаусс, «Theoria motus corporum coelestium» (Гамбург, 1809), излагает тот же метод и добавляет, что пользуется им с 1795 года, то есть с восемнадцатилетнего возраста.

Лежандр был оскорблён. Он написал Гауссу, что о приоритете следует заявлять публикациями, а не воспоминаниями. Гаусс ответил, что метод настолько прост, что многие наверняка додумались до него самостоятельно, и что его собственные вычислительные записи утрачены. В 1820 году Лежандр напечатал публичное обвинение.

Что известно сегодня. Косвенные подтверждения в пользу Гаусса есть: Ольберс в 1816 году свидетельствовал, что Гаусс показывал ему метод в 1802-м, Бессель — аналогично в 1832-м; в 1799 году вышла публикация с гауссовскими расчётами, где метод, по-видимому, применён. Историки (в частности, Стиглер) склоняются к тому, что оба независимы, Гаусс раньше по времени, Лежандр раньше по публикации — та же формула, что для НьютонаИсаак Ньютонанглийский математик, физик и астроном · 1643–1727Создатель анализа, механики и оптики, четверть века проработавший начальником Монетного двора и потративший на алхимию и богословие больше бумаги, чем на физику. и ЛейбницаГотфрид Вильгельм Лейбницнемецкий философ, математик и дипломат · 1646–1716Придумал знаки $d$ и $\int$, которыми мы пишем анализ до сих пор, — и всю жизнь искал язык, на котором спор можно было бы заканчивать словами «посчитаем»..

Есть и третий претендент: Роберт Адрейн в США опубликовал близкие результаты в 1808 году, независимо от обоих.

Чем ответ Гаусса содержательнее

Карл Фридрих Гаусс. Кристиан Альбрехт Йенсен, 1840
Карл Фридрих Гаусс. Кристиан Альбрехт Йенсен, 1840Christian Albrecht Jensen · Public domain

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

Лежандр обосновал метод удобством. Гаусс задал другой вопрос: при каком законе распределения ошибок метод наименьших квадратов оптимален?

Постановка. Пусть измерения $y_i = \theta + \varepsilon_i$, ошибки независимы и одинаково распределены с плотностью $\varphi$. Оценка максимального правдоподобия максимизирует

$$L(\theta) = \prod_{i=1}^{n}\varphi(y_{i}-\theta).$$

Гаусс постулирует то, что все и так делали: наилучшая оценка неизвестной величины по нескольким измерениям есть среднее арифметическое. И спрашивает: какая $\varphi$ это обеспечивает?

Вывод. Логарифмируя и дифференцируя, условие максимума:

$$\sum_{i=1}^{n} \frac{\varphi'(y_{i}-\theta)}{\varphi(y_{i}-\theta)} = 0.$$

Обозначим $g = (\ln\varphi)' = \varphi'/\varphi$. Требуется, чтобы решением было $\theta=\bar y$ при любых данных $y_i$ и любом $n$. Возьмём $n=2$, $y_1 = a$, $y_2 = -a$: тогда $\bar y = 0$ и условие даёт $g(a)+g(-a)=0$, то есть $g$ нечётна. Взяв $n$ произвольным и варьируя данные, получаем функциональное уравнение

$$g(u+v) = g(u)+g(v),$$

решение которого (при разумных условиях регулярности) линейно: $g(u) = cu$. Значит,

$$\frac{\varphi'}{\varphi} = cu \quad\Longrightarrow\quad \ln\varphi = \frac{cu^{2}}{2}+\text{const} \quad\Longrightarrow\quad \varphi(u) = A\,e^{cu^{2}/2}.$$

Для интегрируемости $c<0$; полагая $c=-1/\sigma^{2}$ и нормируя:

$$\boxed{\;\varphi(u) = \frac{1}{\sigma\sqrt{2\pi}}\,e^{-u^{2}/(2\sigma^{2})}\;}$$

Нормальный закон — не гипотеза, а единственное распределение, при котором среднее арифметическое оптимально. Отсюда и название «гауссово».

И обратно: если ошибки нормальны, то максимум правдоподобия

$$\prod_i e^{-(y_i-f_i)^2/2\sigma^2} = \exp\left(-\frac{1}{2\sigma^{2}}\sum_i (y_i-f_i)^{2}\right)$$

достигается при минимуме суммы квадратов. Метод наименьших квадратов и нормальный закон — две стороны одного утверждения.

ЛапласПьер-Симон Лапласфранцузский математик, астроном и физик · 1749–1827Свёл небесную механику в пять томов, а теорию вероятностей — в один, и на полтора столетия задал образ мира, в котором будущее вычисляется из настоящего. в 1810-м замкнул круг с другой стороны (Лаплас): центральная предельная теорема объясняет, почему ошибки измерений в самом деле часто нормальны — они складываются из множества мелких независимых причин.

Церера — с оговоркой

Знаменитый сюжет: 1 января 1801 года Пиацци в Палермо открывает малую планету Цереру, наблюдает её 41 день (дуга всего в 9°), после чего она теряется в лучах Солнца. Найти её заново по такой короткой дуге никто не может. Двадцатичетырёхлетний Гаусс вычисляет орбиту и указывает область неба, сильно отличающуюся от прогнозов остальных. 7 декабря 1801 года Церера найдена фон Цахом, 1 января 1802-го — Ольберсом, вблизи гауссовского положения. Слава мгновенная.

Оговорка, которую стоит проговорить. Часто пишут, что Гаусс «применил метод наименьших квадратов и нашёл Цереру». Историки (Тиц и Уайтхед и другие) относятся к этому осторожнее: основную роль сыграл его собственный новый метод определения орбиты по трём наблюдениям, а наименьшие квадраты, по-видимому, использовались при уточнении по всем данным Пиацци. Сам Гаусс в дневнике связывает метод с возвращением «беглянки», но однозначного технического подтверждения нет.

Продолжение

Замечание для карты

Это точка с двумя равноправными местами: Париж (Лежандр, 1805) и Гёттинген (Гаусс, 1809). Как переписка ПаскаляБлез Паскальфранцузский математик, физик и религиозный философ · 1623–1662Собрал в шестнадцать лет теорему о шестиугольнике, в девятнадцать — счётную машину, в тридцать один — теорию вероятностей, а потом бросил математику ради Бога. и ФермаПьер Фермафранцузский юрист и математик · 1607–1665Советник тулузского парламента, не напечатавший при жизни ни одной математической книги — и успевший заложить теорию чисел, аналитическую геометрию, метод касательных и теорию вероятностей. (Париж — Тулуза), сюжет не рассказывается из одного города. Если схема допускает связь «многие-ко-многим» между событием и местом, это хороший тестовый случай — можно даже добавить третью привязку (Адрейн, США, 1808).

Открыть на карте