691 KiB
Глава 3. Теория оценок¶
3.1 Понятие о точечной оценке¶
Математическая статистика имеет огромное количество разообразнейших применений, но с точки зрения экспериментальной физики (и как следствие студентов, изучающих эту науку) наиболее интересным применением является оценка параметров закономерностей. Пусть есть некоторое явление природы, которое можно описать при помощи модели $M(\theta)$. Здесь $\theta$ - это некоторый набор параметров модели, которые могут принимать различные значения. На этом этапе мы не оговариваем, как именно модель описывает процесс и что мы можем принимать в качестве параметров. Положим теперь, что существует некоторый выделенный набор параметров $\theta_0$, который соответствует некоторому "истинному" состоянию природы. Далее мы будем исходить из того предположения, что при попытке предпринять некоторые измерения, мы будем получать результаты, соответствующие нашей модели именно с этим набором параметров.
Предаставим теперь, что мы провели некоторую серию экспериментов $X = \{X_0, X_1,...X_N\}$, в которых мы тем или иным способом изучаем состояние природы (будем дальше называть результаты этих экспериментов экспериментальной выборкой). Нашей задачей в этой главе будет описание процедуры, при помощи которой можно на основе выборки сделать вывод об истинном состоянии природы $\theta_0$. Важно понимать, что в общем случае, результаты измерений являются случайными величинами, поэтому полученное нами на основании этих данных состояние природы также будет случайной величиной в противовес истинному состоянию природы $\theta_0$, которое вообще говоря, истинно случайной величиной не является. Полученную величину будем называть точечной оценкой состояния природы $\hat{\theta}$ или просто оценкой. Саму процедуру, в процессе которой получена оценка, будем называть оцениванием.
В дальнейшем будем считать, что состояния природы описываются действительным числом или набором действительных чисел. Сама по себе теория этого не требует, но в противном случае довольно сложно сравнивать состояния между собой (требуется определять понятие близости в произвольном пространстве). В этом случае наша процедура оценивания: \begin{equation} \hat{\theta} = f(X) \end{equation} является действительной функцией на пространстве векторов $X$, состоящих из случайных переменных. Такие функции еще называют статистиками. Очевидно, что далеко не людая такая функция будет давать тот результат, которого мы хотим. Поэтому вводятся дополнительные обязательные свойства оценок.
3.1.1 Свойства точечных оценок¶
Состоятельность¶
Естественное пожелание к оценщику, заключается в том, что качество оценки должно зависеть от объема выборки, числа измерений $n$ случайных переменных $X$: чем больше выборка, тем качественней оценка $\hat{\theta}$. Иными словами, мы хотим, чтобы с ростом объема выборки значение оценки приближалось к истинному значению параметра. При использовании сходимости по вероятности оценку $\hat{\theta}$ определяют как состоятельную, если при любых $\varepsilon > 0$ и $\eta > 0$, найдется такое $N$, что вероятность $P \left( \left| \hat{\theta} - \theta \right| > \varepsilon \right) < \eta $ при всех $n > N$. Или другими словами, всегда можно выбрать такое количество измерений, для которого вроятность отклонения оценки от истинного значения не превышает наперед заданное число.
Несмещенность¶
Рассмотрим набор измерений, каждое из которых состоит из $n$ наблюдений $X$, характеризуемый функцией плотности вероятности $P(X | \theta) = P(\hat\theta | \theta)$ при фиксированном $n$ и определим смещение как отклонение среднего по этому набору $\hat{\theta_n}$ от истинного \begin{equation} b = E[\hat{\theta_n}] - \theta \end{equation} Оценка называется несмещенной, если $b = 0$.
Заметим, что смещение не зависит от измеренных величин, но зависит от размера образца, формы оценщика и от истинных (в общем случае неизвестных) свойств ФПВ $f(x)$, включая истинное значение параметра. Если смещение исчезает в пределе $n \to \infty$, говорят об асимптотически несмещенной оценке. Заметим, что из состоятельности оценки не следует несмещенность. Это означает, что даже если $\hat{\theta}$ сходится к истинной величине $\theta$ в единичном эксперименте с большим числом измерений, нельзя утверждать, что среднее $\hat{\theta}$ по бесконечному числу повторений эксперимента с конечным числом измерений $n$ будет сходится к истинному $\theta$. Несмещенные оценки пригодны для комбинирования результатов разных экспериментов. В большинстве практических случаев смещение должно быть мало по сравнению со статистической ошибкой и им пренебрегают.
Эффективность¶
Для сравнения разных методов оценки, очень важным свойством является эффективность. Говоря простым языком, эффективность - это величина, обратная разбросу значений $\hat{\theta}$ при применении к разным наборам данных. Для того, чтобы хорошо разобраться в этом свойстве, надо вспомнить, что оценка, как случайная величина, распределена с плотностью $P(\hat\theta | \theta)$. Вид этого распределения может быть не известен полностью, но знать его свойства – низшие моменты – необходимо. Среднее по нему суть смещение, а дисперсия $\sigma_{\hat\theta}^2 = \int{ (\hat\theta - \theta} ) P(\hat\theta | \theta) d\hat\theta$ суть мера ошибки в определении оценки. Выбирая между различными методами, мы, естественно, хотим, чтобы ошибка параметра была минимальной из всех доступных нам способов его определения для фиксированного эксперимента. Разные методы обладают разной эффективностью и в общем случае при конечной статистике дисперсия распределения оценки никогда не будет равна нулю. Разумеется, встает вопрос о том, можно ли построить оценку с максимальной возможной эффективностью.
3.1.2* Граница Рао-Крамера и информация¶
Утверждение¶
Пусть есть несмещенная оценка $\hat\theta$ параметра $\theta$, тогда всегда выполняется неравенство:
\begin{equation} D(\hat\theta) \geq \frac{1}{I(\theta)}, \end{equation}где $I(\theta)$ - информация Фишера:
\begin{equation} I(\theta) = E_X \left[\left( \frac{\partial \ln L(X,\theta)}{\partial \theta} \right)^2 \right], \end{equation}где в свою очередь $L(X, \theta)$ - функция правдоподобия. $L(X,\theta) = \prod{P(X_i,\theta)}$ - вероятность получить набор данных при фиксированном значении $\theta$.
Доказательство в одномерном случае¶
Введем функцию
\begin{equation} W = \frac{\partial \ln L(X,\theta)} {\partial \theta}. \end{equation}Найдем иатематическое ожидание этой функции:
\begin{equation} E(W) = \int{L(X, \theta) W(X, \theta) dX} = \int{ L \frac{1}{L} \frac{\partial L}{\partial \theta} dX} = \frac{\partial}{\partial \theta} \int{L dX} = 0. \end{equation}Теперь рассмотрим ковариацию $E(\hat\theta W)$:
\begin{equation} E(\hat\theta W) = \int{\hat\theta \frac{1}{L} \frac{\partial L} {\partial {\theta}} L dX} = \frac{\partial}{\partial \theta} \int{\hat\theta L dX} = \frac{\partial}{\partial \theta} E(\hat\theta). \end{equation}Для несмещенных оценок последнее выражение упрощается до вида $E(\hat\theta W) = {\partial \theta}/{\partial \theta} = 1$ Согласно неравенствую Коши — Буняковского, получаем: $\sqrt{D(\hat\theta) D(W)} \geq \left| E(\hat\theta W) \right| = 1$. Отсюда легко получаем желаемое утверждение.
Следствие¶
Максимальная эффективность достигается в том случае, если величины $\hat\theta$ и $W$ являются скоррелированными. Оценка, максимизирующая функцию $L(X,\theta)$ является в общем случае состоятельной, несмещенной, кроме того совпадает с оценкой вида $W(\hat\theta, X) = 0$. То есть является максимально эффективной.
3.2 Интервальные оценки¶
На практике применение точечных оценок сильно затруднено тем, что не известно, на сколько каждая такая оценка точна. Действительно, мы можем спокойно утверждать, что слон весит один килограмм если разброс нашей оценки составляет больше массы слона. Для того, чтобы решить эту проблему есть два пути. Первый путь - это на ряду с точечной оценки указывать меру эффективности этой оценки или ее разброс$\sigma_{\hat\theta}$. Но тут любой внимательный слушатель заметит, что для определения эффективности, вообще говоря, надо знать истинное значение парематра $\theta$, которого мы, разумеется не знаем. Следовательно приходится использовать не эффективность, а оценку этой эффективности, которая сама по себе является случайной величиной. Кроме того, часто случается, что распределение оценки является не симметричным и описать его одним числом не удается.
Более корректным способом является построение интервальной оценки (доверительного интервала). Формально определение интервальной оценки будет отличаться в зависимости от того, какое определение вероятности вы будете использовать.
Частотная интерпретация: интервальной оценкой параметра или группы параметров $\theta$ с уровнем достоверности $\alpha$ называется такая область на пространстве параметров (в одномерном случае - интервал), которая при многократном повторении эксперимента с вероятностью (частотой) $\alpha$ перекрывает истинное значение $\theta$.
Субъективная интерпретация: доверительным интервалом для параметров $\theta$ будем называть такую область в пространстве параметров, в которой интегральная апостериорная вероятность нахождения истинного значения параметра равна $\alpha$.
Для точного описания результата проведения анализа как правило в качестве результата приводят как точечную оценку, так и интервальную оценку с некоторым уровнем достоверности (в английском варианте Confidence Level или C. L.). В некоторых случаях приводят несколько интервальных оценок с разным уровнем достоверности. В случае, когда речь идет об определении верхней или нижней границы какого-то параметра, точечная оценка как правило не имеет смысла и в качестве результата дается только итнервальная оценка.
3.3 Методы построения оценок¶
3.3.1 Метод максимума правдоподобия¶
Определение¶
Введем функцию правдоподобия следующим образом: $L(X,\theta) = \prod{P(X_i,\theta)}$, где $P(X_i, \theta)$ - вероятность получить случайным образом компоненту данных с номером $i$ при истинном значении параметра $\theta$. Будем называеть оценкой максимума правдоподобия такую $\hat\theta$, для которой $L(\hat\theta, X)$ максимально для заданного набора данных. Такая оценка в общем случае является состоятельной и несмещенной. Она имеет очень хорошо понятный физический смысл: $\hat\theta$ - это такое значение параметра, для которого вероятность получить набор данных максимальна. Кроме того, как упоминалось в параграфе 3.1.2, такая оценка кроме всего прочего будет еще и наиболее эффективной из всех возможных (достаточной статистикой).
Ограничения¶
Оценка является эффективной только если область измереия $X$ не зависит от $\theta$.
Интервальная оценка¶
Для построения интервальной оценки воспользуемся субъективной вероятностью и теоремой Байеса: \begin{equation} P(\theta | X) = \frac{L(X|\theta) \pi(\theta)}{\int{LdX}} \end{equation} Здесь $L$ - функция правдоподобия, а $\pi$ - априорная вероятность для параметра $\theta$. Если нет никакой дополнительной информации оп параметре, то мы можем положить $\pi = 1$. Мы получаем, что распределение значения реального параметра повторяет форму функции правдоподобия. Вероятность того, что параметр $\theta$ лежит в диапазоне от $a$ до $b$ составляет \begin{equation} P(a < \theta < b) = \int_b^a{L(X | \theta)}. \end{equation}
Интервальная оценка в асимптотическом случае¶
Согласно центральной предельной теореме, при достаточно большом количестве данных $X$ (и при разумных предположениях о форме распределения для этих данных), функция правдоподобия будет иметь вид нормального распределения. В этом случае центральный доверительный интервал можно вычислить, пользуясь аналитической формулой. Центральный доверительный интервал с уровнем значимости 68% будет соответствовать диапазону между значениями $\theta$, такими, что значение функции правдоподобия в них отличается от максимума на 0.5.
Логарифмическое правдоподобие¶
В реальных задачах часто работают не с самой функцией правдоподобия, а с ее логарифмом $\ln L(X|\theta) = \sum{\ln P(X_i | \theta)}$. Это связано с тем, что при численных манипуляциях проще работать со сложением, чем с умножением. Очевидно, что максимум правдоподобия будет одновременно и максимумом его логарифма и обратно. Для построения интервальных оценок логарифм прадоподобия не годится.
3.3.2 Метод минимума $\chi^2$¶
Разберем частный, но очень часто встречающийся случай, когда распределение измеренной величины - нормальное (Гауссово): \begin{equation} P(X_i | \theta) = \frac{1}{\sqrt {2 \pi} \sigma_i} e^{- \frac{(X_i - \mu_i(\theta))^2}{2 \sigma_i^2}} \end{equation} Здесь $\mu_i(\theta)$ - модельное значение. Опуская нормировочный множитель, не существенный для нахождения максимума можно записать логарифм функции правдоподобия в виде: \begin{equation} \ln L(X_i | \theta) = - \sum{\frac{(X_i - \mu_i(\theta))^2}{2 \sigma_i^2}} \end{equation} Определим величину \begin{equation} \chi^2 = - 2 \ln L(X_i | \theta) = \sum{\frac{(X_i - \mu_i(\theta))^2}{\sigma_i^2}}. \end{equation} Назовем оценкой минимума $\chi^2$ такое значение параметра, при котором эта величина минимально. Очевидно, что эта оценка эквивалента максимуму правдоподобия для нормально расспределенных измерений. Оценка минимума $\chi^2$ будет состоятельной и несмещенной и для других распределений измерений, но только для таких измерений она будет обладать максимальной эффективностью.
Интервальная оценка¶
Также как и в методе максимума правдоподобия, центральный интервал можно получить из аналитической форму функции правдоподобия (для нормально распределенных ошибок, форма правдоподобия будет нормальной всегда, а не только в асимптотике). Из-за множителя 2, мы получаем, что центральный интервал в 68% будет соответствовать отклонению $\chi^2$ на 1. Отклонение на 2 будет соответствовать уже 95% доверительному интервалу.
Проверка качества фита¶
Дополнительное полезоное свойство суммы $\chi^2$ заключается в том, что она позволяет оценить, насколько хорошо данные описываются моделью. В том случае, когда измерения распределены по нормальному закону и независимы между собой, сумма $\chi^2$ оказывается распределена по одноименному распределению. При хорошем соответствии модели и данных величина $\chi^2 / n $, где $n$ - так называемое количество степеней свободы (количество точек минус количество параметров), должна в среднем быть равна 1. Значения существенно большие (2 и выше) свидетельствуют о плохом соответствии или заниженных погрешностях. Значения меньше 0.5 как правило свидетельствуют о завышенных ошибках.
3.3.3 Метод наименьших квадратов¶
В случа, если все ошибки $\sigma_i$ одинаковы, множитель $\frac{1}{\sigma^2}$ можно вынести за скобку. Для нахождения минимума постоянный множитель не важен, поэтому мы можем назвать оценкой наименьших кавдратов такое значение параметра $\hat\theta$, при котором миниальна сумма квадратов: \begin{equation} Q = \sum{(X_i - \mu_i(\theta))^2}. \end{equation} Эта оценка очевидно сохраняет все свойства оценки минимума $\chi^2$, правда только в том случае, если все ошибки действительно одинаковы.
Оценка наименьших квадратов удобна в том случае, когда не известны ошибки отдельных измерений. Она является состоятельной и асимптотически несмещенной, хотя ее эффективность оптимальна только для ограниченного количества случаев. В некоторых случаях, ошибку измерения можно оценить по разбросу данных, используя критерий Пирсона.
Интервальная оценка¶
Оценка методом наименьших квадратов очевидно игнорирует информацию о погрешностях измерений и не позволяет напрямую оценить погрешность результата без дополнительных предположений. Для получения оценки погрешностей надо сделать три предположения:
- Данные описываются предложенной моделью;
- Отклонения данных от модельной кривой независимы между собой (носят статистический характер);
- Статистические ошибки для всех точек равны между собой.
При этих предположениях, можно оценить $\sigma$ для каждой из точек как средне квадратичное отклонение точек от наилучшей модели (той, которая получена минимизацией суммы $Q$). После этого задача получения погрешностей сводится к уже решенной для метода $\chi^2$.
3.3.4 Другие методы¶
Существует довольно большое количество других методов определения параметров зависимости. Как правило они специально заточены под конкретную задачу. В качестве примера можно привести метод моментов. Методо дает состоятельную и несмещенную, но довольно слабо эффективную оценку, но при этом позволяет очень быстро делать оценки для больших объемов данных. Еще один пример - метод квази-оптимальных весов Ф. В. Ткачева, напротив обладает высокой эффективностью и стабильностью, но несколько сложнее в реализации.
3.4 Многопараметрические оценки¶
Однопараметрические оценки очень просты для понимания и реализации, но довольно редко встречаются на практике. Даже при оценке параметров линейно зависимости вида $y = k x + b$ уже существует два параметра: $k$ - наклон прямой и $b$ - смещение. Все перечисленные выше математические методы отлично работают и в многомерном случае, но процесс поиска экстремума функции (максимума в случае метода максимума правдоподобия и минимума в случае методов семейства наименьших квадратов) и интерпретация результатов требуют использования специальных программных пакетов.
3.4.1 Доверительные области в многомерном случае¶
Принцип построения доверительной области в многомерном случае точно такой же, как и для одномерных доверительных интервалов. Требуется найти такую областью пространства параметров $\Omega$, для которой вероятностное содержание для оценки параметра $\hat \theta$ (или самого параметра $\theta$ в засимости от того, какой философии вы придерживаетесь) будет равно некоторой наперед заданной величине $\alpha$: \begin{equation} P(\theta \in \Omega) = \int_\Omega{L(X | \theta)}d\Omega = \alpha. \end{equation}
Реализация на практике этого определения сталкивается с тремя проблемами:
- Взятие многомерного интеграла от произвольной функции - не тривиальная задача. Даже в случае двух параметров, уже требуется некоторый уровень владения вычислительной математикой и компьютерными методами. В случае большего числа параметров, как правило надо использовать специально разработанные для этого пакеты.
- Определить центральный интервал для гипер-области гораздо сложнее, чем сделать это для одномерного отрезка. Единых правил для выбора такой области не существует.
- Даже если удалось получить доверительную область, описать такой объект в общем случае не так просто, так что представление результатов составляет определенную сложность.
Для решения этих проблем, пользуются следующим приемом: согласно центральной предельной теореме, усреднение большого количества одинаково распределенных величин дает нормально распределенную величину. Это же верно и в многомерном случае. В большинстве случаев, мы ожидаем, что функция правдоподобия будет похожа на многомерное нормальное распределение: \begin{equation} L(\theta) = \frac{1}{(2 \pi)^{n/2}\left|\Sigma\right|^{1/2}} e^{-\frac{1}{2} (x - \mu)^T \Sigma^{-1} (x - \mu)} \end{equation} где n - размерность вектора параметров, $\mu$ - вектор наиболее вероятных значений, а $\Sigma$ - ковариационная матрица распределения.
Для многомерного нормального распределения, линии постоянного уровня (то есть поверхности, на которых значение плотности вероятности одинаковые) имеют вид гипер-эллипса, определяемого уравнением $(x - \mu)^T \Sigma^{-1} (x - \mu) = const$. Для любого вероятностного содержания $\alpha$ можно подобрать эллипс, который будет удовлетворять условию на вероятностное содержание. Интерес правда редставляет не эллипс (в случае размерности больше двух, его просто невозможно отобразить), а ковариацонная матрица. Диагональные элементы этой матрицы являются дисперсиями соответствующих параметров (с учетом всех корреляций параметров).
3.4.2 Аналитическая оценка для линейной модели¶
Линейная модель: $mu = kx + b$
\begin{equation} Q = \sum{(X_i - \mu_i(\theta))^2} = \sum{(y_i - kx_i - b)^2}. \end{equation}Найдем минимум:
\begin{equation} \frac{\partial Q}{\partial k} = \sum{-2x_i y_i +4 k x_i^2 + 2x_i b} = 0 \end{equation}\begin{equation} \frac{\partial Q}{\partial b} = \sum{-2y_i +4 k x_i + 2 b} = 0 \end{equation}Решение: \begin{equation} k = \frac{N \sum{x_i y_i} - \sum{x_i}\sum{y_i}} {N \sum{x_i^2} - \left( \sum{x_i} \right)^2} \end{equation}
\begin{equation} b = \frac{\sum{y_i} - k\sum{x_i}}{N} \end{equation}Аналитическое решеие может быть получено в случае любой полиномиальной модели, но не в общем случае.
Практика¶
Установка и подготовка программного обеспечения¶
При работе с интерактивными примерами предполагается использование интерактивной среды Jupyter. Наиболее простым способом подготовки всего, необходимого для работы является установка дистрибутива Anaocnda (предполагается использование Python версии 3). Для работы необходимы пакеты numpy, pandas и matplotlib, которые входят в поставку Anaconda. Опытные пользователи python могут воспользоваться любым другим способом установки.
Для запуска данного "блокнота" в интерактивном режиме необходимо скопировать его в рабочую директорию на свое усмотрение и запустить команду jupyter notebook в этой директории.
Инструменты¶
Подготовим несколько инструментов для представления данных и отрисовки. Для начала импортируем нужные библиотеки.
import numpy as np # Библиотека для численных методов import pandas as pd # Библиотека для табличного представления данных import matplotlib.pyplot as plt # Библиотека для отрисовки графиков # Команда для интерактивной отрисовки графиков %matplotlib inline # Дополнительный паке plotly для интерактивных графиков from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot import plotly.graph_objs as go init_notebook_mode(connected=True)
Генерация и отрисовка данных¶
Генерация данных по нормальному распределениею с заданными средними и стандартными отклонениями¶
def generate(x, y, err): ''' Функция в качестве данных берет три массива: средние значения по координате x, по координате y и ошибки по координате y. Все массивы должны быть одинаковой длины. Возвращает таблицу pandas.DataFrame, содержащую сгенерированные данные. ''' # Проверяем, что среди переданных аргументов нет отрицательных и нулей if (err <= 0).any(): raise Exception # Генерируем нормально распределенные данные generatedY = y + np.random.normal(scale=err) # Создаем объект типа DataFrame с тремя колонками. Объект в фигурных скобках является словарем return pd.DataFrame({'x' : x, 'y' : generatedY, 'err': err})
Отрисовка функций¶
def plot_function(func, a, b, label="", numPoints = 300): ''' Отрисовать функцию func в интервале (a,b) с заданным количеством точек ''' x = np.linspace(a,b, numPoints, endpoint=True) plt.plot(x, func(x), label=label) def plot_data(data, name=""): ''' Отрисовать данные ''' plt.errorbar(data.x, data.y, data.err, label = name, fmt='o')
Параболический апроксиматор¶
Для нахождения оценок методом наименьших квадратов и обобщенным методом наименьших квадратов ($\chi^2$) потребуется находить минимум функции, которую приближенно можно описать параболической зависимостью. Это можно сделать разными способами, но в одномерном случае это проще всего сделать аппроксимируя параболу по трем произвольным точкам (лучше брать этит точки как можно ближе к минимуму) или используя стандартную библиотечную функцию для полиномиального фита.
class ParabolicFunction: def __init__(self, a=None, b=None, c=None): """ Конструктор по-умолчанию """ self.a = a self.b = b self.c = c def from_point(self, x1, x2, x3, y1, y2, y3): """ Конструктор по трем точкам """ self.a = (y3-(x3*(y2-y1)+x2*y1-x1*y2)/(x2-x1))/(x3*(x3-x1-x2)+x1*x2) self.b = (y2-y1)/(x2-x1) - self.a*(x1+x2) self.c = (x2*y1-x1*y2)/(x2-x1)+self.a*x1*x2 return self def from_function(self, x1,x2,x3, func): """ Апроксиматор по трем точкам """ return self.from_point(x1,x2,x3, func(x1), func(x2), func(x3)) def from_delta_center(self, centr, delta, func): """ Конструктор по центральному значениюи и дельте """ return self.from_function(centr-delta, centr, centr+delta, func) def __call__(self, x): return self.a*x*x+self.b*x+self.c def minX(self): """ Метод класса возвращает абсцису минимального значения """ return -self.b/2/self.a def minY(self): """ Метод класса возвращает ординату минимального значения """ return self(self.minX(self)) def __repr__(self): return '{:f} *x^2 + {:f} *x + {:f}\n'.format(self.a,self.b,self.c) # Статическая функция, определяющая минимум неизвестной функции, считая ее параболой def find_min (est, delta, func): return ParabolicFunction().from_delta_center(est, delta, func).minX()
Проверим работу аппроксиматора:
function = lambda x: (x*x + x)*np.exp(-x/100) x = np.linspace(-4,4) plt.plot(x, function(x), label = "parabola") plt.plot(x, ParabolicFunction().from_function(-1, 0, 1, function)(x), label = "parabola-interpolation") plt.legend();
Данные¶
Генерируем и запоминаем данные¶
В качестве модели берем $y = \alpha x^2,$ где $\alpha$ - параметр модели
# Запоминаем нашу модель. Функция берет в качестве аргумента число a и возвращает другую функцию model = lambda a : lambda x: a*x*x # Создаем модель с фиксированным параметром parabola = model(1.1) # Генерируем и запоминаем данные x = np.arange(1, 10) y = parabola(x) err = np.sqrt(y) data = generate(x, y, err) #строим график plot_data(data) plot_function(parabola, 1, 10)
#Выводим данные data
| x | y | err | |
|---|---|---|---|
| 0 | 1 | 1.985047 | 1.048809 |
| 1 | 2 | 3.762057 | 2.097618 |
| 2 | 3 | 13.843041 | 3.146427 |
| 3 | 4 | 21.982251 | 4.195235 |
| 4 | 5 | 38.273192 | 5.244044 |
| 5 | 6 | 41.333730 | 6.292853 |
| 6 | 7 | 56.615843 | 7.341662 |
| 7 | 8 | 83.232328 | 8.390471 |
| 8 | 9 | 93.295287 | 9.439280 |
Метод наименьших квадратов¶
Вводим оценку наименьших квадратов. Функция воздействует на модель и набор данных и возвращает функцию одного переменного ( $\chi^2(\alpha)$ ).
# Фукнкция np.vectorize превращает обычную функцию языка Python в обобщенную функцию библиотеки numpy def sum_of_squares(model, data): def _sum_of_squares(alpha): np_model = np.vectorize(model(alpha)) return np.sum((data.y - np_model(data.x))**2) return np.vectorize(_sum_of_squares)
Строим график для наименьших квадратов
plot_function(sum_of_squares(model, data), 0.5, 1.5, "Sum of squares")
Находим оценку наименьших квадратов
minsq = find_min(1, 0.1, sum_of_squares(model, data)) print(minsq)
1.2128183841151954
Метод $\chi^2$¶
Вводим оценку модифицированных наименьших квадратов ($\chi^2$)
def chi2(model, data): def _chi2(alpha): f = np.vectorize(model(alpha)) return np.sum(((data.y - f(data.x))/ data.err)**2) return np.vectorize(_chi2)
Строим график
plot_function(chi2(model, data), 0.5, 1.5, "chi2")
Находим оценку $\chi^2$
minchi = find_min(1, 0.1, chi2(model, data)) print(minchi)
1.243237812147532
Строим более детальный график
plot_function(chi2(model, data),minchi - 0.1,minchi + 0.1,"chi2")
Утверждение!¶
Сумма вида $\chi^2 = \sum {\frac{y_i - \mu(x_i)}{\sigma_i^2}}$ в достаточно общем случае распределена по распределению $\chi^2$. Доверительный интервал с достоверностью $1\sigma$ (~68%) можно построить, найдя точки, где эта величина откланяется от минимума на 1. Кроме того, при хорошем согласии эксперимента с теорией величина $\frac{\chi^2}{n-k}$ должна быть порядка 1.
sigma = lambda parabolic : 1/np.sqrt(parabolic.a) print(ParabolicFunction().from_delta_center(1, 0.1,chi2(model, data))) errchi = sigma(ParabolicFunction().from_delta_center(1,0.1,chi2(model, data))) print(errchi)
259.090909 *x^2 + -644.223230 *x + 405.581501 0.06212607441973938
Сравнение оценок¶
Сравниваем оценки
print('Оценка методом наименьших квадратов = {:.3f}'.format(minsq)) print('Оценка методом хи-квадрат = {:.3f} +- {:.3f}'.format(minchi,errchi))
Оценка методом наименьших квадратов = 1.213 Оценка методом хи-квадрат = 1.243 +- 0.062
Строим распределения оценок. Заодно проверяем, что средняя ошибка методом $\chi^2$ соответствует разбросу этой оценки
num = 1000; lqs = np.zeros(num) chi2s = np.zeros(num) model = lambda a: lambda x: a*x**2 alpha = 1 parabola = model(alpha) displq = 0; dispchi2 = 0; errchi = 0; for i in range(num): # генерируем данные (попробовать заменить генератор ошибок на 0.1 * y) x = np.arange(1,10) y = parabola(x) err = y*0.3 #np.sqrt(y) data = generate(x, y, err) # считаем оценку методом наименьших квадратов lqval = find_min(1, 0.1, sum_of_squares(model, data)) #добавляем ее в массив, чтобы потом строить гистограмму lqs[i] = lqval # добавляем к счетчику дисперсии displq += (lqval-1)**2 / num #то же самое для хи-квадрат chi2val = find_min(1, 0.1, chi2(model, data)) chi2s[i] = chi2val dispchi2 += (chi2val-1)**2 / num # Считаем среднюю ошибку методом хи-квадрат errchi += sigma(ParabolicFunction().from_delta_center(1,0.1,chi2(model, data))) / num print("Стандартное отклонение оценки методом наименьших квадратов: {}".format(np.sqrt(displq))) print("Стандартное отклонение оценки методом хи-квадрат: {}".format(np.sqrt(dispchi2))) print("Средняя ошибка методом хи-квадрат: {}".format(errchi)) plt.title("Сравнение разных методов оценок"); plt.hist(lqs, label="least squares", bins=30) plt.hist(chi2s, label="chi2", bins = 30) plt.legend();
Стандартное отклонение оценки методом наименьших квадратов: 0.15713530973827325 Стандартное отклонение оценки методом хи-квадрат: 0.09940335113972994 Средняя ошибка методом хи-квадрат: 0.10000000000000182
lqs_plot = go.Histogram(x=lqs, name = 'least squares') chi2s_plot = go.Histogram(x=chi2s, name = 'chi2') iplot([lqs_plot,chi2s_plot])