Обработка результатов учебного эксперимента 4 Рекомендации по выполнению и представлению результатов работы Список литературы

Глава 5 Приложение

5.1 Корреляции

Напомним, если величины $x$ и $y$ независимы, то среднее значение (математическое ожидание) произведения отклонений $\Delta x=x-\overline{x}$ и $\Delta y=y-\overline{y}$ равно нулю:

\overline{\Delta x\cdot\Delta y}=\overline{\Delta x}\cdot\overline{\Delta y}=0.

Если же $x$ и $y$ не являются полностью независимыми, среднее значение произведения их отклонений может быть использовано как количественная мера их зависимости. Наиболее употребительной мерой зависимости двух случайных величин является коэффициент линейной корреляции:

r_{xy}=\frac{\left<\Delta x\cdot\Delta y\right>}{\sigma_{x}\cdot\sigma_{y}}.

(5.1)

Нетрудно проверить (с помощью неравенства Коши–Буняковского), что $-1\leq r\leq 1$ . В частности, для полностью независимых величин коэффициент корреляции равен нулю, $r=0$ , а для линейно зависимых $y=kx+b$ нетрудно получить $r=1$ при $k>0$ и $r=-1$ при $k<0$ . Примеры промежуточных случаев представлены на рис. TODO.

Замечание. Угловой коэффициент прямой в задаче линейной регрессии (3.7) выражается через коэффициент корреляции как

k=r_{xy}\frac{\sigma_{y}}{\sigma_{x}}.

Если коэффициент $r_{xy}$ близок к единице, говорят, что величины коррелируют между собой (от англ. correlate — находиться в связи).

Отсутствие корреляции $\not\Rightarrow$ независимость.

Отметим, что (2.6) — необходимое, но не достаточное условие независимости величин. На рис. TODO приведён пример очевидно зависимых $x$ и $y$ , для которых $r\approx 0$ .

Корреляция $\not\Rightarrow$ причинность.

Ещё одна типичная ошибка — исходя из большого коэффициента корреляции ( $r\to 1$ ) между двумя величинами сделать вывод о функциональной (причинной) связи между $x$ и $y$ . Рассмотрим конкретный пример. Между током и напряжением на некотором резисторе имеет место линейная зависимость $U=IR$ , и коэффициент корреляции $r_{UI}$ действительно равен единице. Однако обратное в общем случае неверно. Например, ток в резисторе коррелирует с его температурой $T$ , $r_{IT}\to 1$ (больше ток — больше тепловыделение по закону Джоуля–Ленца), однако ясно, что нагрев резистора извне не приведёт к повышению тока в нём (скорее наоборот, так как сопротивление металлов с температурой растёт). Ошибка отождествления корреляции и причинности особенно характерна при исследовании сложных многофакторных систем, например, в медицине, социологии и т.п.

5.2 Свойства точечных оценок

Если измеряется одна физическая величина $x$ , то можно поставить задачу по конечному набору данных $\mathbf{x}=\{x_{i}\}$ ( $i=1\ldots n$ ) оценить параметры случайного распределения, которому подчиняется $x$ . В частности, найти среднее значение (математическое ожидание) $\overline{x}$ и дисперсию $\sigma^{2}$ .

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

\overline{x}\approx\left<x\right>=\frac{1}{n}\sum_{i}x_{i}\qquad

(5.2)

и выборочной дисперсии (1.2):

\sigma^{2}\approx s^{2}_{n}=\frac{1}{n}\sum_{i}(x_{i}-\left<x\right>)^{2}.

(5.3)

Оценка параметров должна давать правильное значение хотя бы в пределе большого числа измерений. Если при $n\to\infty$ оценка стремится к истинному значению параметра,

\lim_{n\to\infty}\hat{\theta}(\mathbf{x})\to\overline{\theta},

то её называют состоятельной. Можно показать (см. [5]), что если у распределения, которому подчиняется случайная величина, существуют конечные средние и дисперсия, то оценки (5.2), (5.3) являются состоятельными.

Несмещенные оценки.

Рассмотрим случай малого числа измерений ( $n\gtrsim 1$ ). Тогда даже если оценка состоятельна, она может давать довольно большую ошибку. При фиксированном $n$ функцию оценки $\hat{\theta}(x_{1},\ldots,x_{n})$ можно рассматривать как случайную величину с некоторым распределением, отличающимся от распределения измеряемой величины. Естественно потребовать, чтобы среднее (математическое ожидание) этого распределения совпадало с истинным значением искомого параметра:

\overline{\hat{\theta}(\mathbf{x})}=\overline{\theta}.

В таком случае оценку называют несмещённой.

Нетрудно показать, что выборочное среднее (5.2) является несмещённой оценкой. А вот оценка $s_{n}^{2}$ из (5.3) таким свойством не обладает. Математическое ожидание для величины $s_{n}^{2}$ при фиксированном $n$ оказывается равно $\overline{s_{n}^{2}}=\frac{n-1}{n}\sigma^{2}$ (предлагаем в качестве упражнения проверить данное утверждение самостоятельно). Именно поэтому при малых $n$ для оценки дисперсии рекомендуется использовать формулу (1.4):

\sigma^{2}\approx s^{2}_{n-1}=\frac{1}{n-1}\sum_{i}(x_{i}-\left<x\right>)^{2}.

Эффективность оценки.

Для сравнения разных методов оценки очень важным свойством является их эффективность. На качественном уровне эффективность — величина, обратная разбросу значений $\hat{\theta}(\mathbf{x})$ при применении к разным наборам данных $\mathbf{x}$ . Как обсуждалось выше, оценка $\hat{\theta}(\mathbf{x})$ есть случайная величина, подчиняющаяся некоторому, в общем случае неизвестному, распределению. Среднее $\overline{\hat{\theta}(\mathbf{x})}$ по этому распределению определяет смещение оценки. А его дисперсия $\sigma^{2}\left(\hat{\theta}\right)$ — как раз мера ошибки в определении параметра. Выбирая между различными методами (минимума хи-квадрат, максимального правдоподобия, наименьших квадратов и т. д.), мы, естественно, хотим, чтобы ошибка была минимальной. Разные статистические методы обладают разной эффективностью и в общем случае при конечном $n$ величина $\sigma^{2}\left(\hat{\theta}\right)$ никогда не будет равна нулю.

Теорема, устанавливающая максимальное значение эффективности оценки, рассмотрена в п. 5.3.

5.3 Максимальная эффективность оценки (граница Рао–Крамера)

Максимальная эффективность оценки ограничена теоремой Рао–Крамера.

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

Пусть оценка $\hat{\theta}$ параметра $\theta$ является несмещённой, тогда всегда выполняется неравенство:

\sigma^{2}(\hat{\theta})\geq\frac{1}{I(\theta)},

(5.4)

где

I(\theta)=\overline{\left(\frac{\partial\ln L}{\partial\theta}\right)^{2}}.

(5.5)

Здесь $L(\mathbf{y},\,\theta)$ — введённая в п. 3.1 функция правдоподобия (вероятность получить набор результатов $\mathbf{y}$ при заданном параметре $\theta$ ). Функцию $I(\theta)$ также называют информацией Фишера.

Доказательство в одномерном случае.

Обозначим

U\equiv\frac{\partial\ln L}{\partial\theta}=\frac{1}{L}\frac{\partial L}{% \partial\theta}

и найдём математическое ожидание этой функции:

\overline{U}=\int U\cdot Ld\mathbf{y}=\int\frac{\partial L}{\partial\theta}d% \mathbf{y}=\frac{\partial}{\partial\theta}\int{Ld\mathbf{y}}=0.

Теперь рассмотрим ковариацию параметра $\theta$ и функции $U$ :

\overline{\hat{\theta}\cdot U}=\frac{\partial}{\partial\theta}\int{\hat{\theta% }Ld\mathbf{y}}=\frac{\partial\overline{\hat{\theta}}}{\partial\theta}.

(5.6)

Для несмещенных оценок математическое ожидание оценки параметра равно самому значению параметра: $\overline{\hat{\theta}}=\theta$ , поэтому последнее выражение есть просто единица. Согласно неравенству Коши–Буняковского имеем

\sigma^{2}(\hat{\theta})\cdot\sigma^{2}(U)\geq\left|\overline{\hat{\theta}% \cdot U}\right|=1,

откуда и следует сделанное утверждение.

Следствие.

Максимальная эффективность достигается в том случае, если величины $\hat{\theta}$ и $U$ коррелируют друг с другом. Оценка, максимизирующая функцию $L(\mathbf{y},\theta)$ (метод максимального правдоподобия), является состоятельной, несмещенной, кроме того совпадает с оценкой вида $U(\mathbf{y},\hat{\theta})=0$ , а значит является максимально эффективной.

5.4 Погрешности коэффициентов построения прямой

Проведём подробный вывод для погрешностей коэффициентов наилучшей прямой $\sigma_{k}$ и $\sigma_{b}$ . Воспользуемся общей формулой (2.11) для погрешности косвенных измерений. Считая, что величины $x_{i}$ известны точно, запишем для погрешности углового коэффициента

\sigma_{k}^{2}=\sum\limits_{i}\left(\frac{\partial k}{\partial y_{i}}\right)^{% 2}\sigma_{y_{i}}^{2}.

Продифференцируем (3.13) по $y_{i}$ :

\frac{\partial k}{\partial y_{i}}=\frac{1}{D_{xx}}\frac{\partial}{\partial y_{% i}}\left(\frac{1}{W}\sum w_{i}x_{i}y_{i}-\left\langle x\right\rangle\frac{1}{W% }\sum w_{i}y_{i}\right)=\frac{w_{i}\left(x_{i}-\left\langle x\right\rangle% \right)}{WD_{xx}},

где $D_{xx}=\left<x^{2}\right>-\left<x\right>^{2}$ , $W=\sum_{i}\sigma_{y_{i}}^{-2}$ . Под угловыми скобками здесь понимается выборочное среднее с весами $w_{i}=1/\sigma_{y_{i}}^{2}$ . Тогда

\sigma_{k}^{2}=\frac{1}{W^{2}D_{xx}^{2}}\sum\limits_{i}w_{i}^{2}\left(x_{i}-% \left\langle x\right\rangle\right)^{2}\sigma_{y_{i}}^{2}.

Учитывая, что $w_{i}\sigma_{y_{i}}^{2}=1$ , получим

\sigma_{k}^{2}=\frac{1}{WD_{xx}}.

(5.7)

Аналогично, для погрешности свободного члена имеем

\sigma_{b}^{2}=\sum_{i}\left(\frac{\partial b}{\partial y_{i}}\right)^{2}% \sigma_{y_{i}}^{2},

где

\frac{\partial b}{\partial y_{i}}=\frac{w_{i}}{W}+\frac{\partial k}{\partial y% _{i}}\left\langle x\right\rangle=\frac{w_{i}}{W}\left(1-\frac{x_{i}-\left% \langle x\right\rangle}{\left\langle x^{2}\right\rangle-\left\langle x\right% \rangle^{2}}\left\langle x\right\rangle\right)=\frac{w_{i}}{W}\frac{\left% \langle x^{2}\right\rangle-x_{i}\left\langle x\right\rangle}{D_{xx}}.

Отсюда, пользуясь (5.7), приходим к формуле (3.11):

\sigma_{b}^{2}=\sigma_{k}^{2}\frac{\left\langle\left(\left\langle x^{2}\right% \rangle-x\left\langle x\right\rangle\right)^{2}\right\rangle}{D_{xx}}=\sigma_{% k}^{2}\left\langle x^{2}\right\rangle.

(5.8)

Случай $\sigma_{y}=\mathrm{const}$ .

В частном случае метода наименьших квадратов (п. 3.6.1), формула (5.7) упрощается:

\sigma_{k}^{2}=\frac{\sigma_{y}^{2}}{nD_{xx}},\qquad\sigma_{b}^{2}=\sigma_{k}^% {2}\left<x^{2}\right>.

(5.9)

Здесь величина $\sigma_{y}$ может быть оценена непосредственно из экспериментальных данных:

\sigma_{y}\approx\sqrt{\frac{1}{n-2}\sum_{i}\Delta y_{i}^{2}},

(5.10)

где $n-2$ — число «степеней свободы» для приращений $\Delta y_{i}=y_{i}-(kx_{i}+b)$ ( $n$ точек за вычетом двух связей (3.13)).

Формул (5.9) и (5.10), вообще говоря, достаточно для вычисления погрешности величины $k$ по известным экспериментальным точкам. Однако часто их объединяют в одно упрощённое выражение. Для этого преобразуем (5.10) следующим образом: учитывая, что $\left\langle y\right\rangle=k\left\langle x\right\rangle+b$ , запишем

\Delta y_{i}=y_{i}-kx_{i}-b=\left(y_{i}-\left\langle y\right\rangle\right)-k% \left(x_{i}-\left\langle x\right\rangle\right).

Возведём в квадрат, усредним и воспользуемся выражением для $k$ в форме (3.9):

\left\langle\Delta y^{2}\right\rangle=D_{yy}+k^{2}D_{xx}-2kD_{xy}=D_{yy}-k^{2}% D_{xx}.

Таким образом,

\sigma_{y}=\sqrt{\frac{n}{n-2}\left(D_{yy}-k^{2}D_{xx}\right)},

и с помощью (5.9) получаем формулы (3.10), (3.11):

\boxed{\sigma_{k}=\sqrt{\frac{1}{n-2}\left(\frac{D_{yy}}{D_{xx}}-k^{2}\right)}% ,\qquad\sigma_{b}=\sigma_{k}\sqrt{\left\langle x^{2}\right\rangle}}.

5.5 Многопараметрические оценки

Однопараметрические оценки просты для понимания и реализации, но относительно редко встречаются на практике. Даже при оценке параметров линейной зависимости $y=kx+b$ требуется уже два параметра: наклон $k$ смещение $b$ . Все рассмотренные выше методы нахождения оптимальных параметров работают и в многомерном случае, но поиск экстремума функций (например, максимума функции правдоподобия или минимума суммы квадратов) и интерпретация результатов требуют, как правило, использования численных методов.

5.5.1 Двумерный случай

Остановимся подробнее на построении прямой. Пусть некоторым методом получены точечные оценки для наилучших значений $\hat{k}$ и $\hat{b}$ . Однако самих значений мало — нас интересует область, в которой могут оказаться параметры $k$ , $b$ с некоторой доверительной вероятностью (например, $P=0,68$ ) — двумерная доверительная область.

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

Если бы $k$ и $b$ были независимы, достаточно было бы найти среднеквадратичные отклонения $\sigma_{k}$ и $\sigma_{b}$ , как это сделано в п. 5.4: тогда искомая доверительная область значений параметров на плоскости $(k,b)$ представляла бы собой эллипс, оси которого параллельны координатным (см. рис. 5.1а).

Однако, если взглянуть, к примеру, на рис. 4.2б, иллюстрирующий графический метод построения прямой, можно убедиться, что при варьировании наклона $k$ обязательно меняется и смещение $b$ . То есть параметры $(k,\,b)$ вообще говоря скореллированы. Количественно отклонения параметров будут характеризоваться ковариационной матрицей:

D=\left(\begin{matrix}D_{kk}&D_{kb}\\ D_{bk}&D_{bb}\end{matrix}\right),

где $D_{kk}=\sigma_{k}^{2}$ , $D_{bb}=\sigma_{b}^{2}$ — дисперсии искомых параметров, а

D_{kb}=D_{bk}=\left<(k-\left<k\right>)\cdot(b-\left<b\right>)\right>=\rho_{kb}% \sigma_{k}\sigma_{b}.

Коэффициент $\rho_{kb}$ называют коэффициентом корреляции и он служит показателем «связанности» параметров. Для полностью независимых параметров он равен нулю, а в случае, если параметры нельзя отличить друг от друга — единице.

По известной теореме линейной алгебры, симметричную матрицу можно привести к диагональному виду поворотом координатных осей. Поэтому доверительная область в таком случае будет представлять собой наклонный эллипс (см. рис. 5.1б), а наклон его осей будет определяться коэффициентом корреляции $r_{kb}$ .

Рис. 5.1: Доверительная область значений коэффициентов прямой а) $k$ и $b$ независимы, б) $k$ и $b$ скоррелированы.

5.5.2 Многомерный случай

Принцип построения доверительной области в многомерном случае точно такой же, как и для одномерных доверительных интервалов. Требуется найти такую областью пространства параметров $\Omega$ , для которой вероятностное содержание для оценки параметра $\hat{\theta}$ будет равно некоторой наперед заданной величине $\alpha$ :

P(\theta\in\Omega)=\int\limits_{\Omega}{L(\mathbf{x}|\theta)}d\Omega=\alpha.

(5.11)

Построение многомерной доверительной области на практике сталкивается с тремя проблемами:

1.

Взятие многомерного интеграла от произвольной функции — не тривиальная задача. Даже в случае двух параметров требуется владение методами вычислительной математики. Соответствующие методы реализованы в специализированных программных пакетах.
2.

Определение центрального интервала для многомерной гиперобласти является неоднозначной задачей.
3.

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

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

L(\theta)=\frac{1}{(2\pi)^{n/2}\left|D\right|^{1/2}}e^{-\frac{1}{2}(\mathbf{x}% -\overline{\mathbf{x}})^{T}D^{-1}(\mathbf{x}-\overline{\mathbf{x}})},

(5.12)

где $n$ — размерность вектора параметров, $\overline{\mathbf{x}}$ — вектор наиболее вероятных значений, а $D$ — ковариационная матрица распределения.

Для многомерного нормального распределения, линии постоянного уровня (то есть поверхности, на которых значение плотности вероятности одинаковые) имеют вид гиперэллипса, определяемого уравнением $(\mathbf{x}-\overline{\mathbf{x}})^{T}D^{-1}(\mathbf{x}-\overline{\mathbf{x}})% =\mathrm{const}$ . Для любого вероятностного содержания $\alpha$ можно подобрать эллипс, который будет удовлетворять условию на вероятностное содержание. Интерес, правда, представляет не сам эллипс (в случае размерности больше двух, его просто невозможно наглядно изобразить), а ковариацонная матрица. Диагональные элементы этой матрицы являются дисперсиями соответствующих параметров (с учетом всех корреляций параметров).

5.5.3 Использование пакета scipy для построение оценки

Существует огромное количество программных пакетов для построения численной оценки параметров. Наиболее доступным и широко используемым является пакет scipy на языке Python. Приведем здесь только пример вызова процедуры оптимизации.

Пусть есть экспериментальные данные, представленные в виде трех колонок: $x$ , $y$ и $e r r$ . Требуется построить наилучшую прямую, описывающую эти данные. Код для этого будет выглядеть следующим образом:


from scipy.optimize import curve_fit
function = lambda x, a, b: a*x + b
popt, pcov = curve_fit(function, xdata = x, ydata = y, sigma = err)

После выполнения этого кода, переменная popt содержит массив из двух значений, соответствующих оценке a и b, а переменная pcov содержит ковариационную матрицу для полученных параметров.

Погрешности параметров можно получить как корни из диагональных элементов ковариационной матрицы:


import numpy as np
sigma_a = np.sqrt(pcov[0,0])
sigma_b = np.sqrt(pcov[1,1])

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