Оценка ковариационных матриц

В статистике иногда ковариационная матрица многомерной случайной величины неизвестна, но должна быть оценена . Оценка ковариационных матриц затем решает вопрос о том, как аппроксимировать фактическую ковариационную матрицу на основе выборки из многомерного распределения . Простые случаи, когда наблюдения являются полными, можно рассматривать с помощью выборочной ковариационной матрицы . Выборочная ковариационная матрица (SCM) является несмещенной и эффективной оценкой ковариационной матрицы, если пространство ковариационных матриц рассматривается как внешний выпуклый конус в R ^{p × p} ; однако, измеренная с использованием внутренней геометрии положительно определенных матриц , SCM является смещенной и неэффективной оценкой. ^[1] Кроме того, если случайная величина имеет нормальное распределение , выборочная ковариационная матрица имеет распределение Уишарта , а ее немного по-другому масштабированная версия является оценкой максимального правдоподобия . Случаи, связанные с отсутствующими данными , гетероскедастичностью или автокоррелированными остатками, требуют более глубокого рассмотрения. Другой проблемой является устойчивость к выбросам , к которым матрицы ковариации выборки весьма чувствительны. ^{[2] [3] [4]}

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

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

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

Дана выборка, состоящая из n независимых наблюдений x ₁ ,..., x _n случайного вектора размерности p x ¹ ( столбец ^- вектор размерности p x 1), несмещенная оценка матрицы ковариации ( p x p )

${\displaystyle \operatorname {\Sigma } =\operatorname {E} \left[\left(X-\operatorname {E} [X]\right)\left(X-\operatorname {E} [X]\right)^{\mathrm {T} }\right]}$

это выборочная ковариационная матрица

${\displaystyle \mathbf {Q} ={1 \over {n-1}}\sum _{i=1}^{n}(x_{i}-{\overline {x}})(x_{i}-{\overline {x}})^{\mathrm {T} },}$

где i -е наблюдение p -мерного случайного вектора, а вектор${\displaystyle x_{i}}$

${\displaystyle {\overline {x}}={1 \over {n}}\sum _{i=1}^{n}x_{i}}$

является выборочным средним . Это верно независимо от распределения случайной величины X , при условии, конечно, что существуют теоретические средние и ковариации. Причина появления фактора n  − 1 вместо n по сути та же самая, что и причина появления того же фактора в несмещенных оценках выборочных дисперсий и выборочных ковариаций , которая связана с тем фактом, что среднее неизвестно и заменяется выборочным средним (см. поправку Бесселя ).

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

${\displaystyle \mathbf {Q_{n}} ={1 \over n}\sum _{i=1}^{n}(x_{i}-{\overline {x}})(x_{i}-{\overline {x}})^{\mathrm {T} }.}$

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

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

При оценке кросс-ковариации пары сигналов, которые являются стационарными в широком смысле , пропущенные выборки не обязательно должны быть случайными (например, допустима подвыборка с произвольным фактором). ^{[ необходима ссылка ]}

Случайный вектор X ∈ R ^p ( «вектор-столбец» размера p ×1) имеет многомерное нормальное распределение с невырожденной ковариационной матрицей Σ точно тогда, когда Σ ∈ R p ^{× p является} положительно определенной матрицей , а функция плотности вероятности X равна

${\displaystyle f(x)=(2\pi )^{-{\frac {p}{2}}}\,\det(\Sigma )^{-{\frac {1}{2}}}\exp \left(-{1 \over 2}(x-\mu )^{\mathrm {T} }\Sigma ^{-1}(x-\mu )\right)}$

где μ ∈ R ^{p ×1} — ожидаемое значение X. Ковариационная матрица Σ — это многомерный аналог того, что в одном измерении было бы дисперсией , и

${\displaystyle (2\pi)^{-{\frac {p}{2}}}\det(\Sigma)^{-{\frac {1}{2}}}}$

нормализует плотность так, чтобы она интегрировалась с 1.${\displaystyle f(x)}$

Предположим теперь, что X ₁ , ..., X _n являются независимыми и одинаково распределенными выборками из распределения выше. На основе наблюдаемых значений x ₁ , ..., x _n этой выборки мы хотим оценить Σ.

Функция правдоподобия имеет вид:

${\displaystyle {\mathcal {L}}(\mu ,\Sigma )=(2\pi )^{-{\frac {np}{2}}}\,\prod _{i=1}^{n}\det(\Sigma )^{-{\frac {1}{2}}}\exp \left(-{\frac {1}{2}}(x_{i}-\mu )^{\mathrm {T} }\Sigma ^{-1}(x_{i}-\mu )\right)}$

Достаточно легко показать, что оценка максимального правдоподобия среднего вектора μ представляет собой вектор « выборочного среднего »:

${\displaystyle {\overline {x}}={\frac {x_{1}+\cdots +x_{n}}{n}}.}$

Подробности см. в разделе об оценке в статье о нормальном распределении ; здесь процесс аналогичен.

Поскольку оценка не зависит от Σ, мы можем просто подставить ее вместо μ в функцию правдоподобия , получив${\displaystyle {\bar {x}}}$

${\displaystyle {\mathcal {L}}({\overline {x}},\Sigma )\propto \det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\sum _{i=1}^{n}(x_{i}-{\overline {x}})^{\mathrm {T} }\Sigma ^{-1}(x_{i}-{\overline {x}})\right),}$

а затем искать значение Σ, которое максимизирует правдоподобие данных (на практике проще работать с log  ).${\displaystyle {\mathcal {L}}}$

Теперь мы переходим к первому удивительному шагу: рассмотрим скаляр как след матрицы 1×1. Это позволяет использовать тождество tr( AB ) = tr( BA ) всякий раз, когда A и B — матрицы, имеющие такую ​​форму, что оба произведения существуют. Мы получаем${\displaystyle (x_{i}-{\overline {x}})^{\mathrm {T} }\Сигма ^{-1}(x_{i}-{\overline {x}})}$

${\displaystyle {\begin{aligned}{\mathcal {L}}({\overline {x}},\Sigma )&\propto \det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\sum _{i=1}^{n}\left(\left(x_{i}-{\overline {x}}\right)^{\mathrm {T} }\Sigma ^{-1}\left(x_{i}-{\overline {x}}\right)\right)\right)\\&=\det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\sum _{i=1}^{n}\operatorname {tr} \left(\left(x_{i}-{\overline {x}}\right)\left(x_{i}-{\overline {x}}\right)^{\mathrm {T} }\Sigma ^{-1}\right)\right)\\&=\det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\operatorname {tr} \left(\sum _{i=1}^{n}\left(x_{i}-{\overline {x}}\right)\left(x_{i}-{\overline {x}}\right)^{\mathrm {T} }\Sigma ^{-1}\right)\right)\\&=\det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\operatorname {tr} \left(S\Sigma ^{-1}\right)\right)\end{выровнено}}}$

где

${\displaystyle S=\sum _{i=1}^{n}(x_{i}-{\overline {x}})(x_{i}-{\overline {x}})^{\mathrm {T} }\in \mathbf {R} ^{p\times p}.}$

${\displaystyle S}$иногда называется матрицей рассеяния и является положительно определенной, если существует подмножество данных, состоящее из аффинно независимых наблюдений (что мы и будем предполагать).${\displaystyle p}$

Из спектральной теоремы линейной алгебры следует , что положительно определенная симметричная матрица S имеет единственный положительно определенный симметричный квадратный корень S ^1/2 . Мы можем снова использовать «циклическое свойство» следа, чтобы записать

${\displaystyle \det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \over 2}\operatorname {tr} \left(S^{\frac {1}{2}}\Sigma ^{-1}S^{\frac {1}{2}}\right)\right).}$

Пусть B = S ^1/2 Σ ⁻¹ S ^1/2 . Тогда выражение выше становится

${\displaystyle \det(S)^{-{\frac {n}{2}}}\det(B)^{\frac {n}{2}}\exp \left(-{1 \over 2}\operatorname {tr} (B)\right).}$

Положительно-определенную матрицу B можно диагонализировать, и тогда задача нахождения значения B , максимизирующего

${\displaystyle \det(B)^{\frac {n}{2}}\exp \left(-{1 \over 2}\operatorname {tr} (B)\right)}$

Поскольку след квадратной матрицы равен сумме собственных значений ( «след и собственные значения» ), уравнение сводится к задаче нахождения собственных значений λ ₁ , ..., λ _{p ,} которые максимизируют

${\displaystyle \prod _{i=1}^{n}\lambda _{i}^{\frac {n}{2}}\exp \left(-{\frac {\lambda _{i}}{2}}\right).}$

Это просто задача исчисления, и мы получаем λ _i = n для всех i. Таким образом, предположим, что Q — матрица собственных векторов, тогда

${\displaystyle B=Q(nI_{p})Q^{-1}=nI_{p}}$

т.е. n раз p × p единичная матрица.

Наконец мы получаем

${\displaystyle \Sigma =S^{\frac {1}{2}}B^{-1}S^{\frac {1}{2}}=S^{\frac {1}{2}}\left({\frac {1}{n}}I_{p}\right)S^{\frac {1}{2}}={\frac {S}{n}},}$

т.е., p × p "выборочная ковариационная матрица"

${\displaystyle {S \over n}={1 \over n}\sum _{i=1}^{n}(X_{i}-{\overline {X}})(X_{i}-{\overline {X}})^{\mathrm {T} }}$

является оценщиком максимального правдоподобия "ковариационной матрицы популяции" Σ. В этот момент мы используем заглавную букву X, а не строчную x, потому что мы думаем о ней "как об оценщике, а не как об оценке", т. е. как о чем-то случайном, распределение вероятностей которого мы могли бы получить, зная его. Можно показать, что случайная матрица S имеет распределение Уишарта с n − 1 степенями свободы. ^[5] То есть:

${\displaystyle \sum _{i=1}^{n}(X_{i}-{\overline {X}})(X_{i}-{\overline {X}})^{\mathrm {T} }\sim W_{p}(\Sigma ,n-1).}$

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

${\displaystyle \ln {\mathcal {L}}(\mu ,\Sigma )=\operatorname {constant} -{n \over 2}\ln \det(\Sigma )-{1 \over 2}\operatorname {tr} \left[\Sigma ^{-1}\sum _{i=1}^{n}(x_{i}-\mu )(x_{i}-\mu )^{\mathrm {T} }\right].}$

Дифференциал этого логарифмического правдоподобия равен

${\displaystyle d\ln {\mathcal {L}}(\mu ,\Sigma )=-{\frac {n}{2}}\operatorname {tr} \left[\Sigma ^{-1}\left\{d\Sigma \right\}\right]-{1 \over 2}\operatorname {tr} \left[-\Sigma ^{-1}\{d\Sigma \}\Sigma ^{-1}\sum _{i=1}^{n}(x_{i}-\mu )(x_{i}-\mu )^{\mathrm {T} }-2\Sigma ^{-1}\sum _{i=1}^{n}(x_{i}-\mu )\{d\mu \}^{\mathrm {T} }\right].}$

Она естественным образом распадается на часть, связанную с оценкой среднего значения, и на часть, связанную с оценкой дисперсии. Условие первого порядка для максимума, , выполняется, когда члены, умножающиеся на и тождественно равны нулю. Предполагая, что (оценка максимального правдоподобия) невырожденна, условие первого порядка для оценки вектора среднего значения имеет вид${\displaystyle d\ln {\mathcal {L}}(\mu ,\Sigma )=0}$${\displaystyle d\mu }$${\displaystyle d\Sigma }$${\displaystyle \Sigma }$

${\displaystyle \sum _{i=1}^{n}(x_{i}-\mu )=0,}$

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

${\displaystyle {\widehat {\mu }}={\bar {X}}={1 \over n}\sum _{i=1}^{n}X_{i}.}$

Это позволяет нам упростить

${\displaystyle \sum _{i=1}^{n}(x_{i}-\mu )(x_{i}-\mu )^{\mathrm {T} }=\sum _{i=1}^{n}(x_{i}-{\bar {x}})(x_{i}-{\bar {x}})^{\mathrm {T} }=S}$

как определено выше. Тогда термины, включающие в, могут быть объединены как${\displaystyle d\Sigma }$${\displaystyle d\ln L}$

${\displaystyle -{1 \over 2}\operatorname {tr} \left(\Sigma ^{-1}\left\{d\Sigma \right\}\left[nI_{p}-\Sigma ^{-1}S\right]\right).}$

Условие первого порядка будет выполнено, когда член в квадратных скобках равен (матрично-значному) нулю. Предварительное умножение последнего на и деление на дает${\displaystyle d\ln {\mathcal {L}}(\mu ,\Sigma )=0}$${\displaystyle \Sigma }$${\displaystyle n}$

${\displaystyle {\widehat {\Sigma }}={1 \over n}S,}$

что, конечно, совпадает с каноническим выводом, данным ранее.

Дуайер ^[6] указывает, что разложение на два члена, как показано выше, «необязательно» и выводит оценщик в двух рабочих строках. Обратите внимание, что может быть нетривиально показать, что такой выведенный оценщик является единственным глобальным максимизатором для функции правдоподобия.

Для выборки из n независимых наблюдений x ₁ ,..., x _{n p} -мерной гауссовой случайной величины X с нулевым средним и ковариацией R оценка максимального правдоподобия R определяется как

${\displaystyle {\hat {\mathbf {R} }}={1 \over n}\sum _{i=1}^{n}x_{i}x_{i}^{\mathrm {T} }.}$

Параметр принадлежит множеству положительно-определенных матриц , которое является римановым многообразием , а не векторным пространством , поэтому обычные для векторного пространства понятия ожидания , то есть " ", и смещения оценщика должны быть обобщены на многообразия, чтобы придать смысл проблеме оценки матрицы ковариации. Это можно сделать, определив ожидание многозначной оценки относительно многозначной точки как${\displaystyle R}$${\displaystyle \mathrm {E} [{\hat {\mathbf {R} }}]}$${\displaystyle {\hat {\mathbf {R} }}}$${\displaystyle R}$

${\displaystyle \mathrm {E} _{\mathbf {R} }[{\hat {\mathbf {R} }}]\ {\stackrel {\mathrm {def} }{=}}\ \exp _{\mathbf {R} }\mathrm {E} \left[\exp _{\mathbf {R} }^{-1}{\hat {\mathbf {R} }}\right]}$

где

${\displaystyle \exp _{\mathbf {R} }({\hat {\mathbf {R} }})=\mathbf {R} ^{\frac {1}{2}}\exp \left(\mathbf {R} ^{-{\frac {1}{2}}}{\hat {\mathbf {R} }}\mathbf {R} ^{-{\frac {1}{2}}}\right)\mathbf {R} ^{\frac {1}{2}}}$

${\displaystyle \exp _{\mathbf {R} }^{-1}({\hat {\mathbf {R} }})=\mathbf {R} ^{\frac {1}{2}}\left(\log \mathbf {R} ^{-{\frac {1}{2}}}{\hat {\mathbf {R} }}\mathbf {R} ^{-{\frac {1}{2}}}\right)\mathbf {R} ^{\frac {1}{2}}}$

являются экспоненциальным отображением и обратным экспоненциальным отображением соответственно, «exp» и «log» обозначают обычную матричную экспоненту и матричный логарифм , а E[·] является обычным оператором ожидания, определенным в векторном пространстве, в данном случае касательном пространстве многообразия. ^[1]

Поле собственного вектора смещения оценщика SCM определяется как${\displaystyle {\hat {\mathbf {R} }}}$

${\displaystyle \mathbf {B} ({\hat {\mathbf {R} }})=\exp _{\mathbf {R} }^{-1}\mathrm {E} _{\mathbf {R} }\left[{\hat {\mathbf {R} }}\right]=\mathrm {E} \left[\exp _{\mathbf {R} }^{-1}{\hat {\mathbf {R} }}\right]}$

В этом случае внутреннее смещение оценки определяется как .${\displaystyle \exp _{\mathbf {R} }\mathbf {B} ({\hat {\mathbf {R} }})}$

Для сложных гауссовых случайных величин это векторное поле смещения можно показать ^[1], что оно равно

${\displaystyle \mathbf {B} ({\hat {\mathbf {R} }})=-\beta (p,n)\mathbf {R} }$

где

${\displaystyle \beta (p,n)={\frac {1}{p}}\left(p\log n+p-\psi (n-p+1)+(n-p+1)\psi (n-p+2)+\psi (n+1)-(n+1)\psi (n+2)\right)}$

и ψ(·) — дигамма-функция . Внутреннее смещение матрицы ковариации выборки равно

${\displaystyle \exp _{\mathbf {R} }\mathbf {B} ({\hat {\mathbf {R} }})=e^{-\beta (p,n)}\mathbf {R} }$

и SCM асимптотически несмещен при n → ∞.

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

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

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

Простая версия оценки сжатия матрицы ковариации представлена ​​оценкой сжатия Ледуа-Вольфа. ^{[7] [8] [9] [10]} Рассматривается выпуклая комбинация эмпирической оценки ( ) с некоторой подходящей выбранной целью ( ), например, диагональной матрицей. Затем выбирается параметр смешивания ( ) для максимизации ожидаемой точности сжатой оценки. Это можно сделать с помощью перекрестной проверки или с помощью аналитической оценки интенсивности сжатия. Можно показать, что полученная регуляризованная оценка ( ) превосходит оценку максимального правдоподобия для небольших выборок. Для больших выборок интенсивность сжатия уменьшится до нуля, следовательно, в этом случае оценка сжатия будет идентична эмпирической оценке. Помимо повышенной эффективности оценка сжатия имеет дополнительное преимущество, заключающееся в том, что она всегда положительно определена и хорошо обусловлена.${\displaystyle A}$${\displaystyle B}$${\displaystyle \delta }$${\displaystyle \delta A+(1-\delta )B}$

Были предложены различные цели сокращения:

1. матрица идентичности , масштабированная по средней выборочной дисперсии ;
2. одноиндексная модель ;
3. модель постоянной корреляции, в которой выборочные дисперсии сохраняются, но все парные коэффициенты корреляции предполагаются равными друг другу;
4. двухпараметрическая матрица, где все дисперсии идентичны, а все ковариации идентичны друг другу (хотя и не идентичны дисперсиям);
5. диагональная матрица, содержащая выборочные дисперсии по диагонали и нули во всех остальных местах;
6. матрица идентичности . ^[8]

Оценщик усадки можно обобщить до многоцелевого оценщика усадки, который использует несколько целей одновременно. ^[11] Программное обеспечение для вычисления ковариационного оценщика усадки доступно в R (пакеты corpcor ^[12] и ShrinkCovMat ^[13] ), в Python ( библиотека scikit-learn [1]) и в MATLAB . ^[14]

• Распространение неопределенности
• Выборочное среднее и выборочная ковариация
• Компоненты дисперсии