Uppskattning av kovariansmatriser

I statistik är ibland kovariansmatrisen för en multivariat slumpvariabel inte känd utan måste uppskattas . Uppskattning av kovariansmatriser behandlar sedan frågan om hur man approximerar den faktiska kovariansmatrisen på basis av ett urval från multivariatfördelningen . Enkla fall, där observationerna är fullständiga, kan hanteras med hjälp av provets kovariansmatris . Sampelkovariansmatrisen (SCM) är en opartisk och effektiv estimator ^{av × p} kovariansmatrisen om utrymmet av kovariansmatriser ses som en yttre konvex kon i Rp ; Men mätt med hjälp av den inneboende geometrin hos positiv-definita matriser är SCM en partisk och ineffektiv skattare. Dessutom, om den slumpmässiga variabeln har en normalfördelning , har provets kovariansmatris en Wishart-fördelning och en något annorlunda skalad version av den är den maximala sannolikhetsuppskattningen . Fall som involverar saknade data , heteroskedasticitet eller autokorrelerade residualer kräver djupare överväganden. En annan fråga är robustheten mot extremvärden , för vilka samvariansmatriser är mycket känsliga.

Statistiska analyser av multivariat data involverar ofta explorativa studier av hur variablerna förändras i förhållande till varandra och detta kan följas upp av explicita statistiska modeller som involverar variablernas kovariansmatris. Således spelar uppskattningen av kovariansmatriser direkt från observationsdata två roller:

• att tillhandahålla initiala uppskattningar som kan användas för att studera de inbördes sambanden;
• för att ge provuppskattningar som kan användas för modellkontroll.

Uppskattningar av kovariansmatriser krävs i de inledande stadierna av huvudkomponentanalys och faktoranalys , och är också involverade i versioner av regressionsanalys som behandlar de beroende variablerna i en datamängd, tillsammans med den oberoende variabeln som resultatet av ett slumpmässigt urval .

Uppskattning i ett allmänt sammanhang

Givet ett urval bestående av n oberoende observationer x ₁ ,..., x _n av en p -dimensionell slumpmässig vektor X ∈ R ^{p ×1} (a p ×1 kolumnvektor), en opartisk skattare av ( p × p ) kovariansmatris

$\operatörsnamn {\Sigma } =\operatörsnamn {E} \left[\left(X-\operatörsnamn {E} [X]\höger)\left(X-\operatörsnamn {E} [X]\höger)^{\mathrm {T} }\höger]$

är provets kovariansmatris

$\mathbf {Q} ={1 \över {n-1}}\summa _{ i=1}^{n}(x_{i}-{\overline {x}})(x_{i}-{\overline {x}})^{\mathrm {T} },$

där $x_{i}$ är den i -te observationen av den p -dimensionella slumpmässiga vektorn, och vektorn

${\overline {x}}={1 \över {n}}\summa _{i=1}^{n}x_{i}$

är provmedelvärdet . Detta gäller oavsett fördelningen av stokastisk variabel X , givetvis förutsatt att de teoretiska medelvärdena och kovarianserna existerar. Orsaken till faktorn n − 1 snarare än n är i huvudsak densamma som orsaken till att samma faktor förekommer i opartiska skattningar av urvalsvarianser och urvalskovarianser , vilket relaterar till det faktum att medelvärdet inte är känt och ersätts av urvalet betyda (se Bessels rättelse ).

I de fall fördelningen av den slumpmässiga variabeln X är känd för att vara inom en viss familj av distributioner, kan andra uppskattningar härledas på basis av det antagandet. Ett välkänt exempel är när den slumpmässiga variabeln X är normalfördelad : i det här fallet skiljer sig den maximala sannolikhetsskattaren för kovariansmatrisen något från den opartiska skattningen, och ges av

$\mathbf {Q_{n}} ={1 \över n}\summa _{i=1}^{n}(x_{i}-{\overline {x}})(x_{i}- {\overline {x}})^{\mathrm {T} }.$

En härledning av detta resultat ges nedan. Det är klart att skillnaden mellan den opartiska skattaren och den maximala sannolikhetsuppskattaren minskar för stort n .

I det allmänna fallet ger den opartiska uppskattningen av kovariansmatrisen en acceptabel uppskattning när datavektorerna i den observerade datamängden alla är kompletta: det vill säga att de inte innehåller några saknade element . Ett tillvägagångssätt för att uppskatta kovariansmatrisen är att behandla uppskattningen av varje varians eller parvis kovarians separat och att använda alla observationer för vilka båda variablerna har giltiga värden. Om man antar att de saknade data saknas slumpmässigt resulterar detta i en uppskattning för kovariansmatrisen som är opartisk. För många applikationer kanske denna uppskattning inte är acceptabel eftersom den uppskattade kovariansmatrisen inte garanteras vara positiv semidefinitiv. Detta kan leda till uppskattade korrelationer som har absoluta värden som är större än ett och/eller en icke-inverterbar kovariansmatris.

Vid uppskattning av korskovariansen för ett par signaler som är stationära med bred känsla behöver saknade sampel inte vara slumpmässiga (t.ex. är subsampling med en godtycklig faktor giltig).

Maximal sannolikhetsuppskattning för den multivariata normalfördelningen

Huvudartikel: Multivariat normalfördelning

En slumpvektor X ∈ R ^p (a p ×1 "kolumnvektor") har en multivariat normalfördelning med en icke-singular kovariansmatris Σ exakt om Σ ∈ R ^{p × p} är en positiv-definitiv matris och sannolikhetstäthetsfunktionen för X är

$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)$

där μ ∈ R ^{p ×1} är det förväntade värdet på X . Kovariansmatrisen Σ är den flerdimensionella analogen till vad som i en dimension skulle vara variansen , och

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

normaliserar densiteten $f(x)$ så att den integreras till 1.

Antag nu att X ₁ , ..., X _n är oberoende och identiskt fördelade sampel från fördelningen ovan. Baserat på de observerade värdena x ₁ , ..., x _n i detta prov vill vi uppskatta Σ.

Första stegen

Sannolikhetsfunktionen är:

${\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)$

Det är ganska lätt att visa att den maximala sannolikhetsuppskattningen av medelvektorn μ är vektorn " provmedelvärde" :

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

Se avsnittet om uppskattning i artikeln om normalfördelningen för detaljer; processen här är liknande.

Eftersom uppskattningen ${\bar {x}}$ inte beror på Σ, kan vi bara ersätta μ i sannolikhetsfunktionen ,

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

och sök sedan värdet på Σ som maximerar sannolikheten för data (i praktiken är det lättare att arbeta med log ${\displaystyle {\mathcal {L}}} )$ .

Spåret av en 1 × 1 matris

Nu kommer vi till det första överraskande steget: betrakta skalären $(x_{i}-{\overline {x}})^{\ mathrm {T} }\Sigma ^{-1}(x_{i}-{\overline {x}})$ som spåret av en 1×1 matris. Detta gör det möjligt att använda identiteten tr( AB ) = tr( BA ) närhelst A och B är matriser så formade att båda produkterna existerar. Vi får

${\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 \över 2}\summa _{i=1}^{n}\operatörsnamn {tr} \left(\left(x_{ i}-{\overline {x}}\right)\left(x_{i}-{\overline {x}}\right)^{\mathrm {T}}\Sigma ^{-1}\right)\ höger)\\&=\det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \över 2}\operatörsnamn {tr} \left(\summa _{ 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 }\operatörsnamn {tr} \left(S\Sigma ^{-1}\right)\right)\end{aligned}}$

var

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

$S$ kallas ibland för spridningsmatrisen , och är positiv definitivt om det finns en delmängd av data som består av $p$ affint oberoende observationer (vilket vi kommer att anta).

Med hjälp av spektralsatsen

Det följer av spektralsatsen för linjär algebra att en positiv-definitiv symmetrisk matris S har en unik positiv-definitiv symmetrisk kvadratrot S ^1/2 . Vi kan återigen använda spårets "cykliska egenskap" för att skriva

$\det(\Sigma )^{-{\frac {n}{2}}}\exp \left(-{1 \över 2}\operatörsnamn {tr} \left(S^{\frac {1 }{2}}\Sigma ^{-1}S^{\frac {1}{2}}\right)\right).$

Låt B = S ^1/2 Σ ⁻¹ S ^1/2 . Då blir uttrycket ovan

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

Den positiva-definita matrisen B kan diagonaliseras, och sedan problemet med att hitta värdet på B som maximerar

$\det(B)^{\frac {n}{2}}\exp \left(-{1 \over 2}\ operatörsnamn {tr} (B)\right)$

Eftersom spåret av en kvadratmatris är lika med summan av egenvärden ( "spår och egenvärden" ), reduceras ekvationen till problemet med att hitta egenvärdena λ ₁ , ..., λ _p som maximerar

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

Detta är bara ett kalkylproblem och vi får λ _i = n för alla i. Antag alltså att Q är matrisen av egenvektorer

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

dvs n gånger p × p identitetsmatrisen.

Avslutande steg

Äntligen får vi

$\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}},$

dvs p × p "sampel kovariansmatris"

${S \over n}={1 \over n}\summa _{i=1}^ {n}(X_{i}-{\överlinje {X}})(X_{i}-{\överlinje {X}})^{\mathrm {T} }$

är estimatorn för maximal sannolikhet för "populationskovariansmatrisen" Σ. Vid denna tidpunkt använder vi ett stort X snarare än ett litet x eftersom vi tänker på det "som en estimator snarare än som en uppskattning", dvs. som något slumpmässigt vars sannolikhetsfördelning vi skulle kunna tjäna på att veta. Slumpmatrisen S kan visas ha en Wishart-fördelning med n − 1 frihetsgrader. Det är:

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

Alternativ härledning

En alternativ härledning av skattaren för maximal sannolikhet kan utföras via matriskalkylformler (se även differential för en determinant och differential för den inversa matrisen ). Den verifierar också det tidigare nämnda faktumet om den maximala sannolikhetsuppskattningen av medelvärdet. Skriv om sannolikheten i loggformuläret med spårningstricket:

$\ln {\mathcal {L}}(\mu ,\Sigma )=\operatörsnamn {const} -{n \over 2}\ln \det(\Sigma )-{1 \over 2}\operatörsnamn { tr} \left[\Sigma ^{-1}\sum _{i=1}^{n}(x_{i}-\mu )(x_{i}-\mu )^{\mathrm {T} } \höger].$

Skillnaden mellan denna log-sannolikhet är

$d\ln {\mathcal {L}}(\mu ,\Sigma )=-{\frac {n}{2}}\operatörsnamn {tr} \left[\Sigma ^{-1}\left\ {d\Sigma \right\}\right]-{1 \över 2}\operatörsnamn {tr} \left[-\Sigma ^{-1}\{d\Sigma \}\Sigma ^{-1}\summa _{i=1}^{n}(x_{i}-\mu )(x_{i}-\mu )^{\mathrm {T} }-2\Sigma ^{-1}\summa _{i =1}^{n}(x_{i}-\mu )\{d\mu \}^{\mathrm {T} }\right].$

Det delas naturligtvis upp i den del som är relaterad till uppskattningen av medelvärdet och till den del som är relaterad till uppskattningen av variansen. Första ordningens villkor för maximum, ${\displaystyle d\ln {\mathcal {L}}(\mu ,\Sigma )=0} ,$ är uppfyllt när termerna multiplicerar ${ \displaystyle d\mu }$ och $d\Sigma$ är identiskt noll. Om man antar (den maximala sannolikhetsuppskattningen av) $\Sigma$ är icke-singular, är första ordningens villkor för skattningen av medelvektorn

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

vilket leder till den maximala sannolikhetsskattaren

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

Detta låter oss förenkla

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

som definierats ovan. Då kan termerna som involverar $d\Sigma$ i $d\ln L$ kombineras som

$-{1 \over 2}\operatörsnamn {tr} \left(\Sigma ^{-1}\left\{d\Sigma \right\}\left[nI_{p}-\Sigma ^{-1 }S\höger]\höger).$

Första ordningens villkor $d\ln {\mathcal {L}}(\mu ,\Sigma )=0$ kommer att gälla när termen inom hakparentesen är (matris- värderad) noll. Att förmultiplicera det senare med $\Sigma$ och dividera med $n$ ger

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

vilket givetvis sammanfaller med den kanoniska härledningen som tidigare givits.

Dwyer påpekar att nedbrytning i två termer som visas ovan är "onödig" och härleder estimatorn i två arbetslinjer. Observera att det kanske inte är trivialt att visa att en sådan härledd estimator är den unika globala maximeraren för sannolikhetsfunktion.

Inre kovariansmatrisuppskattning

Inneboende förväntningar

Givet ett urval av n oberoende observationer x ₁ ,..., x _n av en p -dimensionell noll-medelvärde Gaussisk stokastisk variabel X med kovarians R , ges den maximala sannolikhetsestimatorn för R av

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

Parametern $R$ tillhör uppsättningen av positiv-definita matriser , som är en Riemann-manifold , inte ett vektorrum , därav de vanliga vektor-rymdsföreställningarna om förväntan , dvs. " $\ mathrm {E} [{\hat {\mathbf {R} }}]$ ", och estimatorbias måste generaliseras till mångfalder för att förstå problemet med kovariansmatrisuppskattning. Detta kan göras genom att definiera förväntan av en grenrörsvärderad estimator ${\hat {\mathbf {R} }}$ med avseende på den grenrörsvärdade punkten $R$ som

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

var

$\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}}$

$\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}}$

är den exponentiella kartan respektive den inversa exponentiella kartan, "exp" och "log" betecknar den ordinarie matrisexponential- och matrislogaritmen , och E[·] är den ordinarie förväntningsoperatorn definierad på ett vektorrum, i detta fall tangentrymden för grenröret.

Bias av provets kovariansmatris

Det inre förspänningsvektorfältet för SCM-estimatorn ${\hat {\mathbf {R} }}$ definieras som

$\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]$

Den inneboende estimatorbiasen ges sedan av $\exp _{\mathbf {R} }\mathbf {B} ({\hat {\mathbf {R} }})$ .

För komplexa Gaussiska slumpvariabler kan detta förspänningsvektorfält visas vara lika

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

var

$\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)\höger)$

och ψ(·) är digammafunktionen . Den inneboende biasen för provets kovariansmatris är lika med

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

och SCM är asymptotiskt opartisk som n → ∞.

På liknande sätt beror den inneboende ineffektiviteten hos provets kovariansmatris på den Riemannska krökningen av utrymmet för positiv-definita matriser.

Uppskattning av krympning

Om urvalsstorleken n är liten och antalet betraktade variabler p är stort, är ovanstående empiriska skattare av kovarians och korrelation mycket instabila. Specifikt är det möjligt att tillhandahålla estimatorer som avsevärt förbättrar den maximala sannolikhetsuppskattningen i termer av medelkvadratfel. Dessutom, för n < p (antalet observationer är mindre än antalet slumpvariabler) blir den empiriska uppskattningen av kovariansmatrisen singular , dvs den kan inte inverteras för att beräkna precisionsmatrisen .

Som ett alternativ har många metoder föreslagits för att förbättra uppskattningen av kovariansmatrisen. Alla dessa tillvägagångssätt bygger på begreppet krympning. Detta är implicit i Bayesianska metoder och i straffade metoder för maximal sannolikhet och explicit i krympningsmetoden av Stein-typ .

En enkel version av en krympningsuppskattare av kovariansmatrisen representeras av Ledoit-Wolfs krympningsuppskattning. Man betraktar en konvex kombination av den empiriska estimatorn ( $A$ ) med något lämpligt valt mål ( $B$ ), t.ex. den diagonala matrisen. Därefter väljs blandningsparametern ( ${\displaystyle \delta } ) för att maximera den förväntade noggrannheten för den krympta estimatorn.$ Detta kan göras genom korsvalidering eller genom att använda en analytisk uppskattning av krympningsintensiteten. Den resulterande regulariserade estimatorn ( ${\displaystyle \delta A+(1-\delta )B})$ kan visas överträffa den maximala sannolikhetsestimatorn för små prover. För stora prover kommer krympningsintensiteten att minska till noll, så i detta fall kommer krympningsuppskattaren att vara identisk med den empiriska skattaren. Förutom ökad effektivitet har krympningsuppskattningen den ytterligare fördelen att den alltid är positiv och väl konditionerad.

Olika krympmål har föreslagits:

1. identitetsmatrisen , skalad med den genomsnittliga urvalsvariansen ;
2. enkelindexmodellen ; _
3. konstantkorrelationsmodellen, där urvalsvarianserna bevaras, men alla parvisa korrelationskoefficienter antas vara lika med varandra;
4. matrisen med två parametrar, där alla varianser är identiska och alla kovarianser är identiska med varandra (även om de inte är identiska med varianserna);
5. den diagonala matrisen innehåller provvarianser på diagonalen och nollor överallt annars;
6. identitetsmatrisen . _

Krympningsuppskattaren kan generaliseras till en flermålskrympningsuppskattare som använder flera mål samtidigt. Mjukvara för att beräkna en kovarianskrympningsuppskattare finns i R (paket corpcor och ShrinkCovMat ), i Python ( scikit-learn- bibliotek [1] ) och i MATLAB .

Se även

• Spridning av osäkerhet
• Urvalsmedelvärde och urvalskovarians
• Varianskomponenter