Привести оценку матрицы ковариации ошибок

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

2.3.1 Постановка задачи получения оптимальных линейных оценок.

Итак, будем полагать, что задача оценивания (2.1.20), (2.1.21) решается в условиях, когда оцениваемый вектор x и ошибки измерения v считаются случайными. Количественная характеристика качества оценивания может быть введена с помощью скалярной функции

− ~

L(x x ( y)) , устанавливающей определенный штраф за отличие оценки от истинного значения

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

квадратичная функция потерь, имеющая вид

L(x − x ( y)) = (x − x ( y))

− x ( y))

= ∑(xi − xi ( y))

i=1

Введем связанный с этой функцией потерь критерий в виде

( y))

(2.3.1)

J = M x,y {(x − x

(x − x ( y))}= M x,y {SpP},

где

( y))

(2.3.2)

P = M x,y (x

− x ( y))(x − x

матрица ковариаций ошибок ε( y) = x − x ( y) .

Сформулируем задачу оценивания вектора

по измерениям (2.1.21)

следующим образом:

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

xˆ( y) = arg min M	n	~	2
x, y ∑	(x − x ( y))		.
~		i
x	( y)	i =1

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

обеспечивающая его минимум – оптимальной в среднеквадратическом смысле оценкой. Такую оценку называют еще оценкой с минимальной дисперсией (minimum variance estimate).

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

Введем составной случайный вектор z = (xт , y т )т , включающий вектор оцениваемых

параметров и вектор измерений. Из выражения для критерия (2.3.1) следует, что для его

вычисления в общем случае требуется располагать совместной ф.п.р.в. fx,y (x, y) . Обычно при

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

совместной ф.п.р.в. fx,v (x, v) . При наличии fx,v (x, v) плотность fx,y (x, y) в принципе может

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

достаточные основания для введения той или иной ф.п.р.в. fx,v (x, v) , а во-вторых, даже при наличии fx,v (x, v) нахождение fx,y (x, y) не является тривиальной задачей.

С целью упрощения задачи будем считать, что имеется информация лишь о первых двух

моментах вектора z , представленная в виде его математического ожидания z = (x т	, y т )т и
матрицы ковариаций
	z		P x			P xy
P	=		xy		т		,	(2.3.3)
		)	P	y
		(P
а при нахождении оценки используется линейное относительно измерений выражение вида
~		= x + K ( y	− y) .		(2.3.4)
x ( y)

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

~ =

M y x ( y) x .

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

Подставляя (2.3.4) в (2.3.2), можем записать

J = M x,y {(x − x + K ( y − y))т (x − x + K ( y − y))}=

{ } (2.3.5)

= M x,y Sp(x − x + K ( y − y))(x − x + K ( y − y))т .

Из этого выражения вытекает, что для вычисления математического ожидания в (2.3.5) теперь уже достаточно знать только z = (x т , y т )т и P z . Введенные ограничения позволяют сформулировать следующую задачу.

Располагая математическим ожиданием и матрицей ковариации для вектора z = (x т , y т )т ,

найти несмещенную оценку вида (2.3.4), обеспечивающую минимум среднеквадратического критерия (2.3.5). Поскольку предполагается линейная зависимость оценки от измерений, то речь

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

2.3.2 Решение задачи нахождения оптимальных линейных оценок

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

Для того, чтобы линейная оценка (2.3.4) при решении задачи (2.1.20), (2.1.21)

обеспечивала минимум критерия (2.3.5) необходимо и достаточно, чтобы матрица	K opt ,
используемая при вычислении этой оценки, удовлетворяла уравнению
K opt P y = P xy .	(2.3.6)

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

(2.3.4) можно использовать выражение x ( y) = Ky .

Необходимость. Представим произвольную матрицу

K в виде

(2.3.7)

K = K opt

+δK ,

где K opt — матрица, удовлетворяющая (2.3.6),

а δ

— малый скалярный параметр. Подставляя

(2.3.7) в (2.3.5) и раскрывая скобки, получим

opt

J = M x,y {Sp(x −

+ δK ) y)(x −

+ δK ) y)

= Sp[P x

− K opt P yx − P xy (K opt )т

+ K opt (K opt )т ]−

+ P

xy ~ т

~ y

opt

)

− K

opt

y ~ т

)+ δ

~ y ~ т

− δSp(KP

− KP

SpKP K

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

= Sp[P

− 2P

opt

)

+ K

opt

)

] −

− 2δSp(P

xy ~ т

− K

opt

y ~ т

) + δ

SpKP

Поскольку матрица K opt по предположению минимизирует выбранный критерий, то должно выполняться условие

~			~	~			~

	dJ
			=	2Sp(P xy K т − K opt P y K	т )= 2Sp(P xy − K opt P y )K т = 0 .
	dδ
		δ=0
						~	K opt должна
Очевидно, что для выполнения этого условия при	любой	матрице
K ,
удовлетворять уравнению (2.3.6).
Достаточность
Предположим теперь, что K opt удовлетворяет (2.3.6).	Покажем, что оценка	xˆ( y) = K opt y
минимизирует критерий (2.3.5).
Подставляя эту оценку в (2.3.5), запишем

J opt =( M x,y {Sp(x − K opt y)(x − K opt y)т}=).

= Sp P x − 2P xy (K opt )т + K opt P y (K opt )т

Для произвольной матрицы K , заданной в виде (2.3.7), получим

~	~	~ ~	т .
J	= J opt − 2δSp((P xy − K opt P y )K	т )+ δ2SpKP y K

Поскольку по предположению второе слагаемое обращается в ноль, а третье — в силу

неотрицательности P y неотрицательно, то

~
J opt ≤ J ,
что и завершает доказательство достаточности.
Если матрица P y не вырождена, то тогда из (2.3.6) следует
K opt = P xy (P y )−1 ,	(2.3.8)
и таким образом
xˆ( y) = Pxy (P y )−1 y.	(2.3.9)
Сняв требование о центрированном характере векторов x и	y , нетрудно убедиться в том, что

приведенное доказательство полностью сохранится и на этот случай, если в качестве оптимальной

оценки принять выражение (2.3.4), в котором K = K opt , т.е.
xˆ( y) = x + K opt ( y − y),	(2.3.10)
где K — матрица соответствующей размерности.
Используя (2.3.10), нетрудно получить выражение для матрицы ковариаций
Popt = M x, y{(x − xˆ( y))(x − xˆ( y))T } =
= Px + K opt P y (K opt )T − K opt P yx − Pxy (K opt )T .
Преобразуя это выражение с использованием (2.3.8), получаем
Popt = P x − P xy (P y )-1 P yx .	(2.3.11)

Таким образом, алгоритм вычисления оптимальных в среднеквадратическом смысле линейных оценок и соответствующей им матрицы ковариаций при решении задачи оценивания (2.1.20), (2.1.21) задается соотношениями (2.3.10), (2.3.11).

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

оценок	полностью определяется первыми двумя моментами для составного вектора
z = (xт	, y т )т и не зависит от вида ф.п.р.в. fx,y (x, y) . Иными словами при произвольном виде

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

2.3.3 Решение линейной задачи

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

z = (x т , y т )т . Как отмечалось выше, при решении задач обработки навигационной информации считаются известными статистические свойства оцениваемого вектора и ошибок измерения. В

нелинейном случае даже при известном виде совместной ф.п.р.в. fx,v (x, v) возникает проблема

нахождения первых двух моментов вектора z = (x т , y т )т . Эта проблема легко решается для линейной задачи (2.1.10), (2.1.11), если заданы соответствующие характеристики для вектора

(x т , v т )т .

B , где матрица B определяет взаимную корреляцию x и v , то

	P	x		P	xy
P z =				=
(P xy )т	P y			.(2.3.12)
				P x		P x H т + B
=			HPx + Bт
			HPx H т + HB + BтH т + R

Принимая во внимание	тот факт, что в этом случае P xy = P x H т	+ B ,
P y = HP x H т + HB + B тH т + R ,	нетрудно конкретизировать выражения (2.3.10), (2.3.11)	(см.

задачу 2.3.1).

В частности, в случае, когда x и v некоррелированы, т.е. B = 0 , соотношения для оптимальной оценки и матрицы ковариаций ее ошибок могут быть представлены как

K opt = (P x H т )(HP x H т + R)−1 ,

Popt = P x − P x H т (HP x H т + R)−1 HP x .

Учитывая справедливость (1.6.27), (1.6.28) можем также записать

Popt = (P x )−1 + H тR−1H −1 ,

K opt = Popt H тR −1 .

(2.3.13)

(2.3.14)

(2.3.15)

(2.3.16)

Интересным представляется сопоставление оптимального в среднеквадратическом смысле

линейного алгоритма с алгоритмами, рассмотренными в предыдущем разделе, т.е. с алгоритмами метода наименьших квадратов и его модификациями. Нетрудно заметить, что при B = 0 выражение (2.3.10) при подстановке в него выражения (2.3.16) и выражение (2.3.15) совпадают с выражениями (2.2.28), (2.2.29), и таким образом полученный алгоритм совпадет с

алгоритмом ММНК, если в его критерии принять D = (P x )−1 , Q = R −1 . Причины такого совпадения обсуждаются в подразделе 2.5.4. Сопоставление между собой различных вариантов

МНК было проведено в подразделе 2.2.4. Поскольку при сделанных предположениях линейная оптимальная оценка совпадает с ММНК, полученные в этом разделе выводы относительно взаимосвязи ММНК с ОМКН и МНК в полном объеме распространяются на оптимальные в среднеквадратическом смысле линейные оценки.

В частности, если считать выполненным условие

(P x )−1 << H тR−1H	(2.3.17)
и, кроме того, принять R = r 2 E , т.е. полагать, что ошибки представляют	собой

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

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

Проиллюстрируем сказанное на примере оценивания коэффициентов полинома первой степени по измерениям (2.1.3), считая, что ошибки измерения являются некоррелированными между собой центрированными случайными величинами с одинаковыми дисперсиями, равными r 2 , а

оцениваемые коэффициенты представляют собой центрированные ( x = 0 ), некоррелированные между собой и с ошибками измерений случайные величины с матрицей ковариаций

1 . Пример реализации измерений типа (2.1.3) приведен на рис. 2.3.1.

σ 12

Используя соотношения (2.3.15), (2.3.16), легко получить следующие выражения для оценок и матрицы ковариаций

σ2

r 2

opt

xˆ

∑ti

r 2

i=1

σ02

r 2

Popt =

∑ti

σ2

r 2

i=1

−1

∑ti

∑yi

i=1

∑ti2

∑ti yi

σ2

r 2

i=1

−1

∑ti

i=1

∑ti2

r 2

i=1

50
40			измеренные значения

30
20
10
0
-10		линейный тренд

-20
-30	1	2	3	4	5	6	7	8	9	10
0
				номер измерения

Рис.2.3.1 Пример реализации ошибок измерений, содержащих линейный тренд.

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

J мнк (x) = ∑( yi − x0 − x1ti )2 ,

(2.3.18)

i=1

и таким образом можно записать

−1

xˆ

мнк

∑ti

∑yi

мнк

∑ti

i=1

= r

i=1

мнк

∑ti

∑ti yi

∑ti

i=1

В ОМНК при диагональной матрице Q в (2.3.18) появятся сомножители qi , и соответственно выражения для оценки и матрицы ковариаций преобразуются к виду

m∑qi

xˆ омнк = mi=1

∑qi ti

i=1

m	−1		m
∑qiti			∑qi yi
i=1			i=1	,
m			m
∑qi ti2		∑qi ti yi
i=1		i=1

m∑qi

Pомнк = mi=1

∑qi ti

i=1

m	−1
∑qi ti
i=1	.
m
∑qi ti2
i=1

Нетрудно заметить, что при qi =	1		, i =	,	эти выражения совпадут с выражениями для
1.m
r 2

		1		0

МНК, а, положив дополнительно D		2		и учитывая, что x = 0 , легко убедиться, что
=		σ0
		0		1

			σ2
					1

оценки ММНК совпадут с оптимальными оценками.

	1		m		1	1	m
Условие (2.3.17) в этом примере сводится к неравенствам		<<		,		<<		∑ti2 , при
σ2	r 2	σ2	r 2
	0				1			j=1

выполнении которых все рассматриваемые оценки будут между собой совпадать.

2.3.4 Свойство ортогональности оптимальных линейных оценок

Докажем следующее весьма важное свойство линейной оптимальной в среднеквадратическом смысле оценки xˆ( y) . В целях простоты сделаем это для случая x = 0 .

Для того, чтобы линейная оценка xˆ( y) = Ky была оптимальной, необходимо и достаточно

выполнение условия
M{(x − xˆ( y)) y т} = M{εy т} = 0,	(2.3.19)

означающего, что ошибка оценки ε = x − xˆ( y) не коррелирована с вектором измерений y .

Когда два вектора некоррелированы, то, как отмечалось в 1.3.5, говорят также об их ортогональности. В связи с этим сформулированное свойство линейной оптимальной оценки называют также свойством или условием ортогональности [1.9].

Необходимость. Пусть оценка оптимальна,	т. е. матрица K = K opt = P xy (P y )−1 ,	определяется
в соответствии с выражением (2.3.8). Отсюда с очевидностью следует, что
M{(x − xˆ( y)) y т} = M{xy т − K opt yy т} = 0.
Достаточность.
Пусть равенство (2.3.19) выполнено. Тогда можем записать
M{(x − Ky) y т} = 0 ,
т.е. P xy = KP y , откуда для матрицы	K	получаем выражение (2.3.8), подтверждающее
оптимальность оценки, удовлетворяющей уравнению (2.3.19).
Приведенное утверждение будет также справедливо, если вместо (2.3.19) записать
M{(x − xˆ( y))xˆ т ( y)} = M{εxˆ т ( y)} = 0,	(2.3.20)
или
M{(x − xˆ( y)(Ly)т} = M{ε(Ly)т} = 0,	(2.3.21)
где L — произвольная матрица размерности n ×m .
В частности, при L = L j = (0,…0,1,0000)т	из (2.3.21) вытекает, что
j
M{εyiт} = 0,	i =		.	(2.3.22)
1.m

Равенства (2.3.19), (2.3.20), (2.3.21) означают тот факт, что ошибка оценки ε = x − xˆ( y)

ортогональна вектору оценки xˆ( y) , вектору измерений y , произвольной линейной

комбинации его компонент и каждой компоненте по отдельности.

Соседние файлы в предмете [НЕСОРТИРОВАННОЕ]

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

Мы предполагаем, что наблюдения независимы и одинаково распределены (iid).

2.6.1. Эмпирическая ковариация

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

Эмпирическая ковариационная матрица выборки может быть вычислена с использованием empirical_covariance функции пакета или путем подгонки EmpiricalCovariance объекта к выборке данных с помощью EmpiricalCovariance.fit метода. Будьте осторожны, чтобы результаты зависели от центрирования данных, поэтому можно использовать assume_centered параметр точно. Точнее, если assume_centered=False, то предполагается, что тестовый набор имеет тот же средний вектор, что и обучающий набор. В противном случае оба должны быть центрированы пользователем и assume_centered=True использоваться.

Примеры

См. В разделе Оценка ковариации усадки: LedoitWolf против OAS и максимальное правдоподобие пример того, как подогнать EmpiricalCovariance объект к данным.

2.6.2. Уменьшение ковариации

2.6.2.1. Базовая усадка

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

В scikit-learn это преобразование (с определяемым пользователем коэффициентом усадки) можно напрямую применить к предварительно вычисленной ковариации с помощью shrunk_covariance метода. Кроме того, сжатая оценка ковариации может быть адаптирована к данным с помощью ShrunkCovariance объекта и его ShrunkCovariance.fit метода. Опять же, результаты зависят от того, центрированы ли данные, поэтому можно использовать assume_centered параметр точно.

Математически это сжатие состоит в уменьшении отношения между наименьшим и наибольшим собственными значениями эмпирической ковариационной матрицы. Это можно сделать, просто сдвинув каждое собственное значение в соответствии с заданным смещением, что эквивалентно нахождению оценщика максимального правдоподобия со штрафом l2 ковариационной матрицы. На практике усадка сводится к простому выпуклому преобразованию: $Sigma_{rm shrunk} = (1-alpha)hat{Sigma} + alphafrac{{rm Tr}hat{Sigma}}{p}rm Id$.

Выбирая величину усадки, $alpha$ сводится к установлению компромисса смещения/дисперсии и обсуждается ниже.

Примеры

См. В разделе Оценка ковариации усадки: LedoitWolf против OAS и максимальное правдоподобие пример того, как подогнать ShrunkCovariance объект к данным.

2.6.2.2. Усадка Ледуа-Вольфа

В своей статье 1 2004 г. О. Ледуа и М. Вольф предлагают формулу для расчета оптимального коэффициента усадки.α который минимизирует среднеквадратичную ошибку между оценочной и реальной ковариационной матрицей.

Оценщик Ледойта-Вольфа ковариационной матрицы может быть вычислен на выборке с ledoit_wolf функцией sklearn.covariance пакета, или ее можно получить иным способом, подгоняя LedoitWolf объект к той же выборке.

Примечание Случай, когда ковариационная матрица популяций изотропна

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

Примеры

См. В разделе Оценка ковариации усадки: LedoitWolf против OAS и максимальное правдоподобие пример того, как подогнать LedoitWolf объект к данным и визуализировать характеристики оценщика Ледуа-Вольфа с точки зрения правдоподобия.

Рекомендации

О. Ледойт и М. Вольф, «Хорошо обусловленная оценка для многомерных ковариационных матриц», Журнал многомерного анализа, том 88, выпуск 2, февраль 2004 г., страницы 365-411.

2.6.2.3. Приблизительная усадка Oracle

Предполагая, что данные распределены по Гауссу, Chen et al. 2 выведена формула, направленная на выбор коэффициента усадки, который дает меньшую среднеквадратичную ошибку, чем та, которая дается формулой Ледуа и Вольфа. Результирующая оценка известна как оценка ковариации Oracle Shrinkage Approximating.

Оценщик OAS ковариационной матрицы может быть вычислен на выборке с помощью oas функции sklearn.covariance пакета, или он может быть получен иным образом путем подгонки OAS объекта к той же выборке.

Компромисс смещения и дисперсии при установке усадки: сравнение вариантов оценок Ледуа-Вольфа и OAS

Рекомендации

Чен и др., «Алгоритмы сжатия для оценки ковариации MMSE», IEEE Trans. на знак. Proc., Volume 58, Issue 10, October 2010.

Примеры

См. В разделе Оценка ковариации усадки: LedoitWolf против OAS и максимальное правдоподобие пример того, как подогнать OAS объект к данным.
См. Оценку Ледуа-Вольфа и OAS для визуализации разницы среднеквадратических ошибок между a LedoitWolf и OAS оценкой ковариации.

2.6.3. Редкая обратная ковариация

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

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

Оценщик GraphicalLasso использует l1 штраф для обеспечения разреженности на высокоточной матрице: тем выше его alpha параметр, тем более разреженной прецизионную матрицу. Соответствующий GraphicalLassoCV объект использует перекрестную проверку для автоматической установки alpha параметра.

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

Примечание Восстановление структуры

Восстановление графической структуры из корреляций в данных — непростая задача. Если вы заинтересованы в таком восстановлении, имейте в виду, что:

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

Математическая формулировка следующая:
$$hat{K} = mathrm{argmin}_K big( mathrm{tr} S K — mathrm{log} mathrm{det} K + alpha |K|_1 big)$$

Где $K$ — матрица точности, которую необходимо оценить, и $S$ — это выборочная ковариационная матрица. $|K|_1$ — сумма модулей недиагональных коэффициентов $K$. Для решения этой проблемы используется алгоритм GLasso из статьи Friedman 2008 Biostatistics. Это тот же алгоритм, что и в glasso пакете R.

Примеры

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

Рекомендации

Фридман и др., «Оценка разреженной обратной ковариации с помощью графического лассо» , Biostatistics 9, стр. 432, 2008 г.

2.6.4. Оценка робастной ковариации

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

В sklearn.covariance пакете реализована надежная оценка ковариации, определитель минимальной ковариации 3 .

2.6.4.1. Определитель минимальной ковариации

Оценщик, определяющий минимальную ковариацию, — это надежный оценщик ковариации набора данных, представленный PJ Rousseeuw в 3 . Идея состоит в том, чтобы найти заданную пропорцию (h) «хороших» наблюдений, которые не являются выбросами, и вычислить их эмпирическую ковариационную матрицу. Эта эмпирическая ковариационная матрица затем масштабируется, чтобы компенсировать выполненный выбор наблюдений («этап согласованности»). Вычислив оценку определителя минимальной ковариации, можно присвоить веса наблюдениям в соответствии с их расстоянием Махаланобиса, что приведет к повторной оценке ковариационной матрицы набора данных («этап повторного взвешивания»).

Rousseeuw и Van Driessen 4 разработали алгоритм FastMCD для вычисления определителя минимальной ковариации. Этот алгоритм используется в scikit-learn при подгонке объекта MCD к данным. Алгоритм FastMCD также одновременно вычисляет надежную оценку местоположения набора данных.

Сырые оценки могут быть доступны как raw_location_ и raw_covariance_ атрибутами из MinCovDet надежного объекта ковариационной оценивани.

Рекомендации

3 PJ Rousseeuw. Наименьшая медиана квадратов регрессии. Дж. Ам Стат Асс, 79: 871, 1984.
4 Быстрый алгоритм для определения определителя минимальной ковариации, 1999, Американская статистическая ассоциация и Американское общество качества, TECHNOMETRICS.

Примеры

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

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

Отказ от указанных двух предпосылок означает, что ковариационная матрица вектора случайных ошибок (таблица, в которой записаны все ковариации между (varepsilon_{i}) и (varepsilon_{j}), см. параграф 3.3) больше не является диагональной матрицей с одинаковыми числами на главной диагонали, как это было в первоначальной классической модели. Теперь ковариационная матрица вектора случайных ошибок Ω — это произвольная ковариационная матрица (разумеется, так как это не совсем любая матрица, а именно ковариационная матрица, то по своим свойствам она является симметричной и положительно определенной).

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

Проанализируем, к каким последствиям приводит отказ от предпосылок №№4-5.

Во-первых, полученная обычным методом наименьших квадратов оценка ({widehat{beta} = left( {X^{‘}X} right)^{- 1}}X’y) остается несмещенной (это свойство мы доказывали, опираясь как раз лишь на первые три предпосылки).

Во-вторых, МНК-оценки хоть и остаются несмещенными, но больше не являются эффективными.

В-третьих, если мы оцениваем ковариационную матрицу вектора оценок коэффициентов (которая нужна для тестирования всевозможных гипотез), то оценка (widehat{V}{{(widehat{beta})} = left( {X^{‘}X} right)^{- 1}}S^{2}) смещена и больше не является корректной.

Чтобы убедиться в этом, посчитаем ковариационную матрицу от (widehat{beta}) в условиях обобщенной модели (при этом мы используем свойства ковариационной матрицы, перечисленные в параграфе 3.3):

(V{left( widehat{beta} right) = V}{leftlbrack {{({X^{‘}X})}^{- 1}X’y} rightrbrack = V}{leftlbrack {{({X^{‘}X})}^{- 1}X'{({mathit{Xbeta} + varepsilon})}} rightrbrack =})

({}V{leftlbrack {{({{({X^{‘}X})}^{- 1}X’})}varepsilon} rightrbrack = left( {X^{‘}X} right)^{- 1}}X^{‘}Vlbrackvarepsilonrbrack{left( {left( {X^{‘}X} right)^{- 1}X^{‘}} right)^{‘} = {({X^{‘}X})}^{- 1}}X^{‘}Omega{left( {left( {X^{‘}X} right)^{- 1}X^{‘}} right)^{‘} = {({X^{‘}X})}^{- 1}}X^{‘}Omega X{({X^{‘}X})}^{- 1})

Так выглядит ковариационная матрица вектора МНК-оценок в обобщенной модели. Ясно, что она не может быть корректно оценена стандартной оценкой (left( {X^{‘}X} right)^{- 1}S^{2}). Следовательно, прежней формулой пользоваться нельзя: если мы будем использовать стандартные ошибки, рассчитанные по обычной формуле (предполагая выполнение предпосылок классической линейной модели), то получим некорректные стандартные ошибки, что может привести нас к неверным выводам по поводу значимости или незначимости тех или иных регрессоров.

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

Поэтому для получения эффективных оценок обычный МНК не подойдет, и придется воспользоваться альтернативным методом — обобщенным МНК (ОМНК, generalized least squares, GLS). Формулу для расчета оценок коэффициентов при помощи ОМНК позволяет получить специальная теорема.

Теорема Айткена

Если

модель линейна по параметрам и правильно специфицирована
({y = {mathit{Xbeta} + varepsilon}},)
матрица Х — детерминированная матрица, имеющая максимальный ранг k,
(E{(varepsilon) = overrightarrow{0}}),
(V{(varepsilon) = Omega}) — произвольная положительно определенная и симметричная матрица,

то оценка вектора коэффициентов модели ({{widehat{beta}}^{} = {({X’Omega^{- 1}X})}^{- 1}}X’Omega^{- 1}y) является:

несмещенной
и эффективной, то есть имеет наименьшую ковариационную матрицу в классе всех несмещенных и линейных по y оценок.

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

Докажем эту теорему.

Из линейной алгебры известно: если матрица (Omega) симметрична и положительно определена, то существует такая матрица P, что

(P’bullet{P = Omega^{- 1}}) ⇔ ({P^{‘} = Omega^{- 1}}bullet P^{- 1})

А раз такое представление возможно, то воспользуемся им для замены переменных. От вектора значений зависимой переменной (y), перейдем к вектору (Pbullet y), обозначив его как вектор ({y^{} = P}bullet y). Аналогичным образом введем матрицу ({X^{} = P}bullet X) и вектор ошибок ({varepsilon^{} = P}bulletvarepsilon).

Вернемся к исходной модели, параметры которой нас и интересуют:

({y = X}{beta + varepsilon})

Умножим левую и правую части равенства на матрицу (P):

({mathit{Py} = mathit{PX}}{beta + mathit{Pvarepsilon}})

С учетом новых обозначений это равенство можно записать так:

({y^{} = X^{}}{beta + varepsilon^{}})

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

Для этого вычислим математическое ожидание нового вектора ошибок:

(E{left( varepsilon^{} right) = E}{left( mathit{Pvarepsilon} right) = P}bullet E{(varepsilon) = P}bullet{overrightarrow{0} = overrightarrow{0}})

Теперь вычислим ковариационную матрицу вектора (varepsilon^{}):

(V{left( varepsilon^{} right) = V}{left( {Pbulletvarepsilon} right) = P}bullet V(varepsilon)bullet{P^{‘} = P}bulletOmegabullet{P^{‘} = P}bulletOmegabulletOmega^{- 1}bullet{P^{- 1} = I_{n}})

Здесь (I_{n}) обозначает единичную матрицу размера n на n.

Следовательно, для модели со звездочками выполняются все предпосылки теоремы Гаусса — Маркова. Поэтому получить несмещенную и эффективную оценку вектора коэффициентов можно, применив к этой измененной модели обычный МНК:

({{widehat{beta}}^{} = left( {{X^{}}^{‘}X^{}} right)^{- 1}}{X^{}}^{‘}y^{})

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

({{widehat{beta}}^{} = left( {{X^{}}^{‘}X^{}} right)^{- 1}}{X^{}}^{‘}{y^{} = left( {{(mathit{PX})}^{‘}mathit{PX}} right)^{- 1}}left( mathit{PX} right)^{‘}{mathit{Py} = left( {X’P’mathit{PX}} right)^{- 1}}X^{‘}P^{‘}{mathit{Py} =}left( {X’Omega^{- 1}P^{- 1}mathit{PX}} right)^{- 1}X’Omega^{- 1}P^{- 1}{mathit{Py} = {({X’Omega^{- 1}X})}^{- 1}}X’Omega^{- 1}y)

Что и требовалось доказать.

Взвешенный МНК, который мы обсуждали ранее, — это частный вариант обобщенного МНК (для случая, когда только предпосылка №4 нарушена, а предпосылка №5 сохраняется).

Как и при использовании взвешенного МНК в ситуации применения ОМНК коэффициент R-квадрат не обязан лежать между нулем и единицей и не может быть интерпретирован стандартным образом.

Слабая сторона ОМНК состоит в том, что для его реализации нужно знать не только матрицу регрессоров X с вектором значений зависимой переменной y, но и ковариационную матрицу вектора случайных ошибок (Omega). На практике, однако, эта матрица почти никогда не известна. Поэтому в прикладных исследованиях практически всегда вместо ОМНК используется, так называемый, доступный ОМНК (его ещё называют практически реализуемый ОМНК, feasible GLS). Идея доступного ОМНК состоит в том, что следует сначала оценить матрицу (Omega) (традиционно обозначим её оценку (widehat{Omega})), а уже затем получить оценку вектора коэффициентов модели, заменив в формуле ОМНК (Omega) на (widehat{Omega}):

({{widehat{beta}}^{} = {({X'{widehat{Omega}}^{- 1}X})}^{- 1}}X'{widehat{Omega}}^{- 1}y.)

Применение этого подхода осложняется тем, что (widehat{Omega}) не может быть оценена непосредственно без дополнительных предпосылок, так как в ней слишком много неизвестных элементов. Действительно, в матрице размер (n) на (n) всего (n^{2}) элементов, и оценить их все, имея всего (n) наблюдений, представляется слишком амбициозной задачей. Даже если воспользоваться тем, что матрица (Omega) является симметричной, в результате чего достаточно оценить только элементы на главной диагонали и над ней, мы все равно столкнемся с необходимостью оценивать (left( {n + 1} right){n/2}) элементов, что всегда больше числа доступных нам наблюдений.

Поэтому процедура доступного ОМНК устроена так:

Делаются некоторые предпосылки по поводу того, как устроена ковариационная матрица вектора случайных ошибок (Omega). На основе этих предпосылок оценивается матрица (widehat{Omega}).
После этого по формуле ({({X'{widehat{Omega}}^{- 1}X})}^{- 1}X'{widehat{Omega}}^{- 1}y) вычисляется вектор оценок коэффициентов модели.

Из сказанного следует, что доступный ОМНК может быть реализован только в ситуации, когда есть разумные основания сформулировать те или иные предпосылки по поводу матрицы (widehat{Omega}). Рассмотрим некоторые примеры таких ситуаций.

Пример 5.4. Автокорреляция и ОМНК-оценка

Рассмотрим линейную модель ({y = X}{beta + varepsilon}), для которой дисперсия случайных ошибок постоянна, однако наблюдается так называемая автокорреляция первого порядка:

(varepsilon_{i} = {{rho ast varepsilon_{i — 1}} + u_{i}})

Здесь (u_{i}) — независимые и одинаково распределенные случайные величины с дисперсией (sigma_{u}^{2}), а (rhoin{({{- 1},1})}) — коэффициент автокорреляции.

(а) Найдите ковариационную матрицу вектора случайных ошибок для представленной модели.

(б) Запишите в явном виде формулу ОМНК-оценки вектора коэффициентов модели, предполагая, что коэффициент (rho) известен.

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

Решение:

(а) Используя условие о постоянстве дисперсии случайной ошибки, то есть условие (mathit{var}{{(varepsilon_{i})} = mathit{var}}{(varepsilon_{i — 1})}), найдем эту дисперсию:

(mathit{var}{{(varepsilon_{i})} = mathit{var}}{({{rho ast varepsilon_{i — 1}} + u_{i}})})

(mathit{var}{{(varepsilon_{i})} = rho^{2}}mathit{var}{{(varepsilon_{i — 1})} + mathit{var}}{(u_{i})})

(mathit{var}{{(varepsilon_{i})} = rho^{2}}mathit{var}{{(varepsilon_{i})} + sigma_{u}^{2}})

(mathit{var}{left( varepsilon_{i} right) = frac{sigma_{u}^{2}}{1 — rho^{2}}}.)

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

(mathit{cov}{left( {varepsilon_{i},varepsilon_{i — 1}} right) = mathit{cov}}{left( {{{rho ast varepsilon_{i — 1}} + u_{i}},varepsilon_{i — 1}} right) = {rho ast mathit{cov}}}{left( {varepsilon_{i — 1},varepsilon_{i — 1}} right) + mathit{cov}}{left( {u_{i},varepsilon_{i — 1}} right) = {rho ast mathit{var}}}{{left( varepsilon_{i — 1} right) + 0} = frac{rho ast sigma_{u}^{2}}{1 — rho^{2}}})

По аналогии легко убедиться, что

(mathit{cov}{left( {varepsilon_{i},varepsilon_{i — k}} right) = frac{rho^{k} ast sigma_{u}^{2}}{1 — rho^{2}} = frac{sigma_{u}^{2} ast rho^{k}}{1 — rho^{2}}}.)

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

({Omega = frac{sigma_{u}^{2}}{1 — rho^{2}}}begin{pmatrix} begin{matrix} 1 & rho rho & 1 end{matrix} & ldots & begin{matrix} rho^{n — 2} & rho^{n — 1} rho^{n — 3} & rho^{n — 2} end{matrix} ldots & ldots & ldots begin{matrix} rho^{n — 2} & rho^{n — 3} rho^{n — 1} & rho^{n — 2} end{matrix} & ldots & begin{matrix} 1 & rho rho & 1 end{matrix} end{pmatrix})

(б) Вектор ОМНК-оценок коэффициентов имеет вид:

({{widehat{beta}}^{} = {({X’Omega^{- 1}X})}^{- 1}}X^{‘}Omega^{- 1}{y =}{left( {X’begin{pmatrix} begin{matrix} 1 & rho rho & 1 end{matrix} & ldots & begin{matrix} rho^{n — 2} & rho^{n — 1} rho^{n — 3} & rho^{n — 2} end{matrix} ldots & ldots & ldots begin{matrix} rho^{n — 2} & rho^{n — 3} rho^{n — 1} & rho^{n — 2} end{matrix} & ldots & begin{matrix} 1 & rho rho & 1 end{matrix} end{pmatrix}^{- 1}X} right)^{- 1} ast}X’begin{pmatrix} begin{matrix} 1 & rho rho & 1 end{matrix} & ldots & begin{matrix} rho^{n — 2} & rho^{n — 1} rho^{n — 3} & rho^{n — 2} end{matrix} ldots & ldots & ldots begin{matrix} rho^{n — 2} & rho^{n — 3} rho^{n — 1} & rho^{n — 2} end{matrix} & ldots & begin{matrix} 1 & rho rho & 1 end{matrix} end{pmatrix}^{- 1}y)

Обратите внимание, что дробь (frac{sigma_{u}^{2}}{1 — rho^{2}}) при расчете представленной оценки сокращается. Поэтому для вычисления оценки знать величину (sigma_{u}^{2}) не нужно.

Примечание: если коэффициент автокорреляции (rho) неизвестен, то его можно легко оценить. Например, для этого можно применить обычный МНК к исходной регрессии, получить вектор остатков и оценить регрессию ({widehat{e}}_{i} = {widehat{rho} ast e_{i — 1}}). Полученной оценки (widehat{rho}) достаточно, чтобы вычислить ОМНК-оценку вектора параметров модели. Тем самым в представленном примере для применения доступного ОМНК достаточно оценить всего один параметр ковариационной матрицы вектора оценок коэффициентов.

Пример 5.5. Гетероскедастичность и ОМНК-оценка

Рассмотрим линейную модель ({y = X}{beta + varepsilon}), для которой выполнены все предпосылки классической линейной модели множественной регрессии за одним исключением: дисперсия случайной ошибки прямо пропорциональна квадрату некоторой известной переменной

(mathit{var}{left( varepsilon_{i} right) = sigma_{i}^{2} = sigma_{0}^{2}}{z_{i}^{2} > 0}.)

(а) Найдите ковариационную матрицу вектора случайных ошибок для представленной модели.

(б) Запишите в явном виде формулу ОМНК-оценки вектора коэффициентов модели.

Решение:

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

(Omega = begin{pmatrix} begin{matrix} {sigma_{0}^{2}z_{1}^{2}} & 0 0 & {sigma_{0}^{2}z_{2}^{2}} end{matrix} & begin{matrix} ldots & 0 ldots & 0 end{matrix} begin{matrix} ldots & ldots 0 & 0 end{matrix} & begin{matrix} ldots & ldots ldots & {sigma_{0}^{2}z_{n}^{2}} end{matrix} end{pmatrix})

(б) Обратите внимание, что при подстановке в общую формулу для ОМНК-оценки величина (sigma_{0}^{2}) сокращается, следовательно, для оценки вектора коэффициентов знать её не нужно:

({{widehat{beta}}^{} = {({X’Omega^{- 1}X})}^{- 1}}X^{‘}Omega^{- 1}{y =}left( {X’begin{pmatrix} begin{matrix} z_{1}^{2} & 0 0 & z_{2}^{2} end{matrix} & begin{matrix} ldots & 0 ldots & 0 end{matrix} begin{matrix} ldots & ldots 0 & 0 end{matrix} & begin{matrix} ldots & ldots ldots & z_{n}^{2} end{matrix} end{pmatrix}^{- 1}X} right)^{- 1}X’begin{pmatrix} begin{matrix} z_{1}^{2} & 0 0 & z_{2}^{2} end{matrix} & begin{matrix} ldots & 0 ldots & 0 end{matrix} begin{matrix} ldots & ldots 0 & 0 end{matrix} & begin{matrix} ldots & ldots ldots & z_{n}^{2} end{matrix} end{pmatrix}^{- 1}y)

* * *

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

Источник

Мы предполагаем, что наблюдения независимы и одинаково распределены (iid).

2.6.1. Эмпирическая ковариация

Примеры

См. В разделе Оценка ковариации усадки: LedoitWolf против OAS и максимальное правдоподобие пример того, как подогнать EmpiricalCovariance объект к данным.