inggeo.fsv.cvut.cz

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
**[[..:04_teorie_chyb|Teorie chyb]]**

**<<** [[04_teorie_chyb:0417_kombinovane_vyrovnani|17. Kombinované vyrovnání]]\\
**>>** [[04_teorie_chyb:0419_metody_reseni_normalnich_rovnic|19. Metody řešení normálních rovnic]]

====== 18. Metody robustního odhadu ======
===== Úvod =====

Klasické statistické postupy zahrnující jak statistické testy, tak vyrovnání metodou nejmenších čtverců jsou velmi důležitou součástí geodézie v oblasti zpracování a analýzy měření. Předpokladem jejich správného fungování je však normální rozdělení chyb, bez splnění této podmínky vytvořený pravděpodobnostní model není správný. Bylo zjištěno, že i malé odchylky od normálního rozdělení pravděpodobnosti mají značný vliv na kvalitu výsledku, to znamená, že i jen několik málo hrubých chyb může znehodnotit jinak kvalitní měření.  

Robustní statistické metody si oproti těm klasickým zachovávají funkčnost v určitém okolí normálního rozdělení, zjednodušeně řečeno neselžou při „mírném" nesplnění požadavku na normální rozdělení chyb, tj. pokud jsou správná měření (vyhovující normálnímu rozdělení) kontaminována odlehlými měřeními. Čím je metoda odolnější oproti vlivu chybných (odlehlých) měření, tím je robustnější. Tyto metody jsou založeny na několika principech, které budou dále ve stručnosti popsány. Je vhodné ještě dodat, že aplikovatelnost většiny metod nepřekračuje využití při výpočtu lineární regrese, přesto je vhodné je v celkovém přehledu uvést. 

Cílem přehledu a následného testování je nalézt metody aplikovatelné na vyrovnání metodou nejmenších čtverců pro případy, kdy počet nadbytečných měření není vzhledem k množství určovaných parametrů příliš velký, příkladem může být vyrovnání volné geodetické sítě v prostoru. Vzhledem k tomu, že při vyrovnání je hledáno „nejlepší" řešení, tj. takové, které poskytuje nejmenší charakteristiky přesnosti, je vhodné využít metody robustní statistiky pro vyhledání hrubých chyb a jejich vyloučení z výpočtu a následně provést vyrovnání metodou nejmenších čtverců. Je známo, že malé množství chybných měření lze odhalit testy odlehlých měření, v případě vyšší kontaminace je vhodné (nutné) použít pro jejich identifikaci dále popsané postupy.

===== Třídy odhadů robustní statistiky =====

Základních metod je několik, v [(Jurečková, J.: //Robustní statistické metody//. Nakladatelství Karolinum, Praha 2001. ISBN 80-246-0259-8.)] lze nalézt nejrozšířenější třídy odhadů reálného parametru a to M-odhady, L-odhady a R-odhady.

M-odhady jsou založeny na metodě maximální věrohodnosti (maximum-likelihood), L-odhady na výpočtu lineárních kombinací pořadových statistik a R-odhady na neparametrických testech. Dále jsou známy také generalizované verze těchto odhadů, případně také S-odhady či $\tau $-odhady. Kromě těchto jsou známy ještě další metody, které lze obtížně zařadit, případně vznikly empiricky, a z těch budou vybrány pouze některé využívané při řešení problémů hraničících s geodézií. 

Pouze některé metody lze prakticky využít pro řešení úloh geodézie, mnohé ostatní lze využít pouze pro určení odhadu polohy rozdělení a tyto nebudou dále zmiňovány. Podle [(Huber, P. J.: //Robust Statistics//. Wiley, New York, 1981.)] umožňují M-odhady dostatečnou flexibilitu a lze je výpočetně snadno zvládnout i v případech řešení generalizovaných lineárních modelů, proto bude následující text zaměřen zejména tímto směrem.

Odvození jsou teoreticky složitá, a proto bude jejich princip podán tak, aby bylo možné zejména pochopit hlavní myšlenku či princip.


===== M-odhady =====

Podle [(Koch, K. R.: //Parameter Estimation and Hypothesis Testing in Linear Models//. Springer Verlag, Berlin Heidelberg New York, 1999. ISBN 3-5406525-74.)] se jedná o řešení minimalizační úlohy. Základem je požadavek na maximální věrohodnost řešení daný výrazem:

$$\overline{{\mathbf x}}{\rm =}{\arg  {\mathop{\sup }_{{\mathbf x}} L\left({\mathbf l}{\rm ,\ }{\mathbf x}\right)\ }\ } ,  $$ 

kde ${\mathbf l}$ je náhodný vektor (měření), $L\left({\mathbf l},\ {\mathbf x}\right)$ věrohodnostní (pravděpodobnostní) funkce, $\overline{{\mathbf x}}$ odhad parametrů (neznámých). Jinak řečeno řešení má být nejpravděpodobnější. Jestliže náhodný vektor pozorování ${\mathbf l}$ má hustotu pravděpodobnosti $f\left({\mathbf x}\right)$, která závisí na fixních a neznámých parametrech ${\mathbf x}$, pak pro věrohodnostní funkci platí:

$$L\left({\mathbf l},\ {\mathbf x}\right)=f\left({\mathbf x}\right) .  $$ 

Lineární model úlohy:

$${\mathbf A}\cdot {\mathbf x}={\mathbf l}+{\mathbf \varepsilon } ,   $$ 

$$A=\left( \begin{array}{c}
{{\mathbf a}}^T_1 \\ 
{{\mathbf a}}^T_2 \\ 
\vdots  \\ 
{{\mathbf a}}^T_n \end{array}
\right) , {\mathbf x}=\left( \begin{array}{c}
x_1 \\ 
x_2 \\ 
\vdots  \\ 
x_u \end{array}
\right) ,$$ 

kde ${\mathbf x}$ je vektor neznámých, ${\mathbf l}$ vektor pozorování (měření), ${\mathbf \varepsilon }$ vektor skutečných chyb, ${\mathbf A}$ matice lineárních (linearizovaných) vztahů mezi měřením a neznámými parametry. Počet měření $n$ je větší než počet neznámých $u$, jde o úlohu s vyrovnáním. Měření jsou nezávislá.

Rovnice pozorování pro měření $l_i$:

$$l_i+{\varepsilon }_i={{\mathbf a}}^T_i\cdot {\mathbf x} .  $$ 

Neznámé ve vektoru ${\mathbf x}$ jsou určovány metodou maximální věrohodnosti, jestliže hustota pravděpodobnosti $p(l_i)$ pozorování $l_i$ je úměrná funkci $f\left(l_i,{\mathbf x}\right)$, tj. platí:

$$p\left(l_i\right)\propto f\left(l_i,{\mathbf x}\right) .  $$ 

Hustota pravděpodobnosti $p({\mathbf l})$ pro nezávislá měření je pak:

$$p\left({\mathbf l}\right)\propto \prod^n_{i=1}{f\left(l_i,{\mathbf x}\right)}=f\left(l_1,{\mathbf x}\right)\cdot f\left(l_2,{\mathbf x}\right)\cdot \dots \cdot f\left(l_n,{\mathbf x}\right) .   $$ 


Pro metodu maximální věrohodnosti je nutný předpoklad znalosti rozdělení pravděpodobnosti. S ohledem na centrální limitní větu se dále předpokládá, že měření mají normální rozdělení $N\left({\mathbf A}\cdot {\mathbf x}{\mathbf ,\ }{\sigma }^2\cdot {\mathbf I}\right)$). Potom pro věrohodnostní funkci platí: 

$$L\left({\mathbf l};{\mathbf x},{\sigma }^2\right)=\frac{1}{{\left(2\pi {\sigma }^2\right)}^{n/2}}\cdot e^{-\frac{{{\mathbf \varepsilon }}^T{\mathbf \varepsilon }}{2{\sigma }^2}}.$$ 

Úkolem je maximalizovat věrohodnost (pravděpodobnost), problém se řeší diferenciací podle neznámých ${\mathbf x}$ a směrodatných odchylek ${\sigma }^2$ a položením získaných výrazů nule. Vzhledem k vlastnostem funkce normálního rozdělení i věrohodnostní funkce lze s výhodou derivovat logaritmus této funkce:

$${\ln  L\left({\mathbf l};{\mathbf x},{\sigma }^2\right)\ }=-\frac{n}{2}{\ln  \left(2\pi \right)\ }-\frac{n}{2}{\sigma }^2-\frac{1}{2{\sigma }^2}\cdot {{\mathbf \varepsilon }}^T{\mathbf \varepsilon }.     $$ 

(Po derivaci podle neznámých ${\mathbf x}$ je výsledkem odhad metody nejmenších čtverců). 

Maximalizuje se tedy výraz

$$\prod^n_{i=1}{f\left(l_i,{\mathbf x}\right)}$$ 

nebo s použitím předchozích logaritmovaných vzorců se minimalizuje

$$\sum^n_{i=1}{\left(-{\ln  f\left(l_i,{\mathbf x}\right)\ }\right)} ,  $$ 

což vede k metodě nejmenších čtverců.

Tato odhadová funkce je dána předpokladem normálního rozdělení měření a musí být nahrazena funkcí vhodnější, protože předpoklad normálního rozdělení není splněn. Výpočet robustních M-odhadů ([(Huber P. J.: //Robust Estimation of a Location Parameter//. Annals of Mathematical Statistics 35(1964), 73-101.)]) je dán minimalizací výrazu:

$$\sum^n_{i=1}{\rho \left(l_i,g_i\left(x\right)\right)}=min., $$ 

kde $\rho $ je vhodná odhadová funkce („score function"), $g_i\left({\mathbf x}\right)$ je funkce neznámých parametů. 

Odhadová funkce není konstantní jako v případě metody nejmenších čtverců. Derivace odhadové funkce (tzv. vlivová funkce):

$$\psi \left(l_i,{\mathbf x}\right)=\frac{\partial }{\partial g_i\left(x\right)}\rho \left(l_i,g_i\left(x\right)\right) .$$ 

Odhad $\overline{{\mathbf x}}$ neznámých parametrů ${\mathbf x}$:

$$\sum^n_{i=1}{\frac{\partial\rho \left(l_i,g_i\left(x\right)\right)}{\partial {\widehat{{\mathbf x}}}_{{\mathbf k}}}}=\sum^n_{i=1}{\psi \left(l_i,\widehat{{\mathbf x}}\right)\cdot {\mathbf \ }}\frac{\partial g_i\left(\widehat{{\mathbf x}}\right)}{\partial {\widehat{{\mathbf x}}}_{{\mathbf k}}}=0\ $$

pro $k\in \left\{1,\dots ,u\right\}$.  

Robustní odhad lze nalézt pomocí funkce $\psi \left(l_i,\overline{{\mathbf x}}\right)$. O robustní odhad se jedná pouze v případě, kdy funkce $\psi \left(l_i,\overline{{\mathbf x}}\right)$ je ohraničená, jelikož je úměrná influenční funkci ([(Hampel, F. R. - Ronchetti, E. M. - Rousseeuw, P. R. - Stahl, W. A.: //Robust Statistics//. Wiley, New York, 1986.)]). Influenční funkce popisuje efekt dalšího pozorování na odhad. Spolu s bodem selhání se jedná o důležité charakteristiky popisující konkrétní robustní odhad. Bod selhání je velmi zjednodušeně řečeno nejmenší podíl pozorování, které po nahrazení libovolnými hodnotami mohou vést k chybným hodnotám odhadu (blíže viz [(#5)]).

Pro běžné výpočty je ještě vhodné zavést pojem normované chyby, tj. podíl chyby a příslušné směrodatné odchylky měření:

$${\widehat{\varepsilon }}_i=\frac{{{\mathbf a}}_i\cdot {\mathbf x}{\mathbf -}l_i}{{\sigma }_i}{\mathbf =}\frac{{\varepsilon }_i}{{\sigma }_i}.  $$ 

Odhadem skutečných chyb jsou opravy $v_i$, normované opravy ${\hat{v}}_i$ jsou definovány:

$${\hat{v}}_i{\mathbf =}\frac{v_i}{{\sigma }_i} .   $$ 

Z toho plyne odhadová funkce:

$$\rho \left(l_i,g_i\left({\mathbf x}\right)\right)=\rho \left({\widehat{\varepsilon }}_i\right) .  $$ 

Protože platí $g_i\left(\overline{{\mathbf x}}\right)={{\mathbf a}}^T_i\cdot \overline{{\mathbf x}}$, lze výpočet odhadu upravit na tvar

$$\sum^n_{i=1}{\psi \left(\frac{{{\mathbf a}}_i\cdot {\mathbf x}{\mathbf -}l_i}{{\sigma }_i}\right)\cdot {\mathbf \ }}\frac{a_{ik}}{{\sigma }_i}=0\ , \sum^n_{i=1}{\psi \left(\frac{v_i}{{\sigma }_i}\right)\cdot \frac{a_{ik}}{{\sigma }_i}}=0.    $$ 

Pro metodu nejmenších čtverců platí:

Hustota pravděpodobnosti normálního rozdělení:

$$\varphi \left(x\right){\rm =}\frac{{\rm 1}}{\sigma \sqrt{{\rm 2}\pi }}\cdot e^{{\rm -}\frac{{\left(x{\rm -}E\left(x\right)\right)}^{{\rm 2}}}{{\rm 2}{\sigma }^{{\rm 2}}}}, \varphi \left(\varepsilon \right){\rm =}\frac{{\rm 1}}{\sigma \sqrt{{\rm 2}\pi }}\cdot e^{{\rm -}\frac{{\varepsilon }^{{\rm 2}}}{{\rm 2}{\sigma }^{{\rm 2}}}},     $$ 

Funkce $\rho $  (bez konstant, po aplikaci logaritmu):

$$\rho \left(\varepsilon \right){\rm =}{\ln  \left(e^{{\rm -}{\frac{{\rm 1}}{{\rm 2}}\left(\frac{\varepsilon }{\sigma }\right)}^{{\rm 2}}}\right)\ }{\rm =-}{\frac{{\rm 1}}{{\rm 2}}\left(\frac{\varepsilon }{\sigma }\right)}^{{\rm 2}} .$$ 

Po derivaci:

$$\psi \left(\frac{\varepsilon }{\sigma }\right){\rm =}\frac{\partial \left({\rm -}{\frac{{\rm 1}}{{\rm 2}}\left(\frac{\varepsilon }{\sigma }\right)}^{{\rm 2}}\right)}{\partial {\rm \ }\left(\frac{\varepsilon }{\sigma }\right)}{\rm =-}\frac{{\rm 1}}{{\rm 2}}\cdot {\rm 2}\cdot \frac{\varepsilon }{\sigma }{\rm =-}\frac{\varepsilon }{\sigma } .     $$ 

Pro minimalizaci musí platit:

$$\sum^n_{i{\rm =1}}{{\rm -}\frac{{\varepsilon }_i}{{\sigma }_i}\cdot \frac{a_{ik}}{{\sigma }_i}}{\mathbf =}\sum^n_{i{\rm =1}}{\frac{{\varepsilon }_i}{{\sigma }_i}\cdot \frac{a_{ik}}{{\sigma }_i}}{\rm =0} .$$ 

Skutečné chyby nejsou známy, budou nahrazeny jejich odhadem - opravami. Řešení normálních rovnic metodou nejmenších čtverců je výpočetně jednoduché. Aby bylo ovšem možné tento výpočet použít, normované opravy ${\hat{v}}_i$ je nutné formálně vynásobit korekčním koeficientem $w_i$. Musí tedy platit:

$$\sum^n_{i{\rm =1}}{\psi \left(\frac{v_i}{{\sigma }_i}\right)\cdot \frac{a_{ik}}{{\sigma }_i}}{\rm =}\sum^n_{i{\rm =1}}{w_i\cdot \frac{v_i}{{\sigma }_i}\cdot \frac{a_{ik}}{{\sigma }_i}}{\mathbf =}0$$ 

a tedy:

$$\psi \left(\frac{v_i}{{\sigma }_i}\right)\cdot \frac{a_{ik}}{{\sigma }_i}{\mathbf =}w_i\cdot \frac{v_i}{{\sigma }_i}\cdot \frac{a_{ik}}{{\sigma }_i} , \Rightarrow  w_i{\mathbf =}\psi \left(\frac{v_i}{{\sigma }_i}\right){\rm /}\frac{v_i}{{\sigma }_i} .    $$ 

Korekční člen $w_i$ představuje určitou váhu měření $l_i$, jejichž velikost je přímo závislá na velikosti normované opravy ${\hat{v}}_i$, tj. $w_i=w_i(v_i,{\sigma }_i)$.

Po zavedení korekčního členu lze provést úpravu:

$$\frac{1}{{\sigma }^2}\cdot \sum^n_{i=1}{{w_i\cdot v}_i\cdot {\mathbf \ }}a_{ik}=0$$

pro $k\in \left\{1,\dots ,u\right\}$.    

Dále je možné definovat diagonální váhovou matici ${\mathbf W}$:

$${\mathbf W}=diag\left(w_1,w_2,\dots ,\ w_n\right)$$ 

a řešit normální rovnice ve tvaru: 

$${{\mathbf A}}^T\cdot {\mathbf W}\cdot {\mathbf A}\cdot {\mathbf x}{\mathbf =}{{\mathbf A}}^T\cdot {\mathbf W}\cdot {\mathbf l} . $$ 

Váhy $w_i$ závisí na opravách $v_i$, tj. na odhadu neznámých ${\mathbf x}$. Odhad tedy musí být určován iterativně, jako první aproximaci lze použít výsledek metody nejmenších čtverců. 

$${\overline{{\mathbf x}}}^{{\mathbf (}0{\mathbf )}}{\mathbf =}{\left({{\mathbf A}}^T\cdot {\mathbf A}\right)}^{{\mathbf -}{\mathbf 1}}\cdot {{\mathbf A}}^T\cdot {\mathbf l}, $$ 

$${\overline{{\mathbf x}}}^{{\rm (}m{\rm +1)}}{\rm =}{\left({{\mathbf A}}^T\cdot {{\mathbf W}}^{{\rm (}m{\rm )}}\cdot {\mathbf A}\right)}^{{\rm -}{\mathbf 1}}\cdot {{\mathbf A}}^T\cdot {{\mathbf W}}^{{\rm (}m{\rm )}}\cdot {\mathbf l} ,    $$ 

$${{\mathbf v}}^{{\mathbf (}m{\mathbf )}}{\mathbf =}{\mathbf A}\cdot {\overline{{\mathbf x}}}^{{\mathbf (}m{\mathbf )}}{\mathbf -}{\mathbf l}. $$ 

V [(#2)] je dokázána konvergence tohoto výpočtu.


==== Huberův M-odhad ====

Jestliže pozorování ${\mathbf l}$ mají normální rozdělení, lze za funkci rozdělení pravděpodobnosti použít (bez konstantních členů):

$$f\left(\varepsilon \right)=e^{-\frac{{\varepsilon }^2}{2}}$$ 

a tedy jako odhadovou funkci použít

$$\rho \left(\varepsilon \right)=\frac{{\varepsilon }^2}{2}.    $$ 

Pak platí 

$$\psi \left(\varepsilon \right)=\frac{\partial\rho \left(\varepsilon \right)}{\partial\varepsilon }=\varepsilon , w_i=\frac{\psi \left(\varepsilon \right)}{\varepsilon }=1,$$ 

a minimalizuje se výraz

$$\frac{{\rm 1}}{{\rm 2}\cdot {\sigma }^2}\sum^n_{i=1}{{\varepsilon }^2_i}.   $$ 

Až na násobný koeficient se jedná o metodu nejmenších čtverců a metodu nerobustní, funkce $\psi \left(\varepsilon \right)=\varepsilon $ neohraničeně roste s rostoucím $\varepsilon $, váha je konstantní $w_i=1$.

Při volbě robustní odhadové funkce vychází Huber ([(#2)]) z normálního rozdělní náhodné veličiny. Jeho řešení je založeno na nahrazení okrajových částí normálního rozdělení pravděpodobnosti Laplaceovým rozdělením (speciální formou exponenciálního rozdělení). Tímto způsobem je dosaženo ohraničení vlivové funkce $\psi \left(\varepsilon \right)$ a dosažení robustnosti odhadu. Toto řešení vede k větší pravděpodobnosti výskytu odlehlých měření na okrajích rozdělení.

{{ :04_teorie_chyb:041801_nerobustni_odhad.jpg }}
;#;
Obr. 1 //(Nerobustní) odhad metodou nejmenších čtverců - tvar funkce// $\psi $
;#;

Pro Huberův odhad platí:

$f\left(\varepsilon \right)=e^{-\frac{{\varepsilon }^2}{2}}$ pro $\left|\varepsilon \right|\le c$ ,      

$f\left(\varepsilon \right)=e^{-c\cdot \left|\varepsilon \right|+c^2/2}$ pro $\left|\varepsilon \right|>c$ ,      

kde $c$ je konstanta závisející na předpokládaném množství kontaminace měřených dat ohlednými hodnotami. Podle [(#2)] se pro 4% odlehlých měření volí $c=1,5$, pro méně než 1% $c=2,0$. Odhadová funkce je tedy dána vztahy:

$\rho \left(\varepsilon \right)=\frac{{\varepsilon }^2}{2}$ pro $\left|\varepsilon \right|\le c$,       

$\rho \left(\varepsilon \right)=c\left|\varepsilon \right|-\frac{c^2}{2}$ pro $\left|\varepsilon \right|>c$ ,  

pro vlivovou funkci $\psi $ platí:

$\psi \left(\varepsilon \right)=\varepsilon$ pro $\left|\varepsilon \right|\le c$,      

$\psi \left(\varepsilon \right)=c\frac{\varepsilon }{\left|\varepsilon \right|}$ pro $\left|\varepsilon \right|>c$ . 

{{ :04_teorie_chyb:041802_huberuv_odhad.jpg }}
;#;
Obr. 2 //Huberův robustní odhad - tvar funkce// $\psi $
;#;     

Funkce $\psi $ je ohraničená, jak je zřejmé z obrázku a tedy se jedná o robustní odhad. Váhy se tedy volí podle předpisů:

$w^{\left(m\right)}_i=1$ pro $\left|v\right|\le c$,      

$w^{\left(m\right)}_i=c\frac{\sigma }{\left|v^{\left(m-1\right)}_i\right|}$ pro $\left|v\right|>c$ .      

Odvození je provedeno pro případ, kdy jsou všechna měření stejně přesná. V případě zavedení vah podle předchozích kapitol je zde uvedená váha $w^{\left(m\right)}_i$ násobný koeficient, kterým se původní váha změní. Rovněž koeficient $c$ hodnotí velikost normované opravy ${\hat{v}}_i$, nikoli opravy s libovolnou variancí.

Huberův M-odhad zde byl prezentován včetně odvození, neboť se jedná o průlomovou práci. Jsou známy další odhadové funkce, jejichž přehled lze je nalézt v dále uvedené literatuře, v českém jazyce je souhrn a potřebná odvození uveden v [(Štroner, M. - Hampacher, M.: //Zpracování a analýza měření v inženýrské geodézii//. 1. vyd. Praha: CTU Publishing House, 2011. 313 s. ISBN 978-80-01-04900-6.)].


<refnotes>
reference-format : []
note-id-format   : []
note-id-base     : text
reference-base   : text
</refnotes>


----
**<<** [[04_teorie_chyb:0417_kombinovane_vyrovnani|17. Kombinované vyrovnání]]\\
**>>** [[04_teorie_chyb:0419_metody_reseni_normalnich_rovnic|19. Metody řešení normálních rovnic]]