Odhad parametrů lineárního modelu metodou nejmenších čtverců (MNČ)

Ilustrační příklad

Naším cílem je najít takovou hodnotu , pro kterou bude platit, že suma čtverců rozdílů mezi jednotlivými měřeními a vypočtenou hodnotou bude minimální. Pro jeden rozdíl můžeme tedy psát

(1) .

(2) .

Cheme-li hodnotu takovou, aby byla funkce M daná vztahem (2) minimální, musíme najít extrém této funkce, kdy bude představovat proměnnou (hodnoty známe z měření jako nějaká konkrétní čísla). Derivací M podle dostáváme

(3) .

(4) ,

a tedy pro hodnotu dostáváme

(5) ,

Pro úplnost uveďme vztahy a postup pro N měření, tedy máme soubor měření , kdy . Poté pro sumu čtverců platí

(6) .

(7) .

Položením vztahu (7) rovno 0 a úpravou dostáváme pro parametr očekávaný výraz pro aritmetický průměr

(8) .

Obecné odvození

MNČ pro měření stejné váhy

(9) ,

kde a popisují fyzikální model měřeného problému, x je vektor parametrů, jejichž hodnoty určujeme zprostředkujícím měřením, a jsou námi hledané parametry, .

(10) ,

, , .

Označíme-li jako vektor parametrů, které budou určeny výpočtem MNČ, poté dosazením do rovnice (10) můžeme napsat vektor odchylek levé a pravé strany této rovnice jako

(11) .

(12) .

Pro maticovou derivaci (12) dle parametrů poté můžeme psát

(13) .

(14) .

Příklad

(15) .

clear variables;

x = linspace( -2, 2, 10 ); y = x;

[ X, Y ] = meshgrid( x, y );

a_0 = 1;

a_1 = 0.5;

a_2 = 0.3;

a_3 = 0.2;

a_4 = -0.3;

Z_nominal = a_0 + a_1*X + a_2*Y + a_3*X.^2 + a_4*Y.^2;

figure

mesh( X, Y, Z_nominal )

xlabel( 'x' )

ylabel( 'y' )

zlabel( 'z = f(x,y)' )

title( 'Nominalní tvar plochy' )

Zavedeme-li do měření z náhodný šum se směrodatnou odchylkou , poté dostáváme vektor měření Z.

Z_mereni = Z_nominal + 0.2*randn( size(Z_nominal) );

figure

mesh( X, Y, Z_mereni )

xlabel( 'x' )

ylabel( 'y' )

zlabel( 'z = f(x,y)' )

title( 'Simulované měření' )

, , ,

kde značí i-té měření z-ové souřadnice, , N je počet všech měření.

MNČ pro měření nestejné váhy (zobecněná metoda nejmenších čtverců)

Předpokládejme nyní, že hodnoty vektoru l jsou různých vah , a jsou vzájemně nezávislé. Metoda nejmenších čtverců je pak definována tak, že minimalizuje sumu čtverců vážených odchylek, tj.

(16) ,

(17)

.

(18) .

Zákon přenášení (šíření) variancí (ZPV)

Základní pojmy

Spojitá jednorozměrná náhodná veličina

Střední hodnota spojité náhodné veličiny X s hustotou pravděpodobnostního rozdělení je definována jako

(19) .

Variance náhodné veličiny X je definována jako

(20) .

Směrodatná odchylka náhodné veličiny X je definována jako

(21) .

Vícerozměrná náhodná veličina

Vektor středních hodnot náhodné veličiny X je obdobně jako v případě (19) definován jako

(22) .

Kovarianční matice (variančně-kovarianční matice) náhodné veličiny X je definována jako

(23) , .

(24)

(25)

značí tzv. kovarianci mezi dvěma náhodnými veličinami a . Kovarianci charakterizuje vzájemnou závislost dvou veličin. Jsou-li tyto veličiny navzájem nezávislé, potom je rovna nule. Jsou-li určitým způsobem závislé, poté nabývá kovariance hodnot od 0 do 1.

(26) .

(27) .

ZPV pro lineární model

Předpokládejme, že máme dva vektory náhodných veličin a , které jsou v následujícím snadném lineárním vztahu

(28) ,

kde A a b jsou konstantní matice a vektor. Jinými slovy každou z veličin můžeme vyjádřit jako lineární kombinaci veličin X. V řadě praktických situacích je toto splněno. Případně můžeme zajistit linearitu např. pomocí Taylorova rozvoje s využitím pouze jeho prvních členů, tzv. linearizaci.

Nechť dále známe kovarianční matici náhodných veličin X. A a b jsou konstanty, tj. nemají kovarianční matici.

Jaká bude v tomto případě kovarianční matice veličin Y? Dosazením do definice kovarianční matice (23) dostáváme

(29) .

Vyjádříme-li střední hodnotu jako

(30) ,

(31)

(32) .

Budou-li veličiny X vzájemně nezávisle, a tedy jejich kovarianční matice bude diagonální, poté pro i-tou veličinu platí

(33) ,

(34) .

ZPV pro nelineární model

(35) ,

kde je vektorová funkce parametrů X, jejichž kovarianční matici známe. Linearizací použítím Taylorova rozvoje a zanedbáním členů vyšších prvního přiblížení (uvážení pouze lineárních členů) dostáváme

(36) ,

kde, značí vektor počátečních parametrů (přibližné hodnoty parametrů, konstanty), dX je vektor přírůstků parametrů, J je tzv. Jakobiho matice prvních derivací. Jestliže jsou počáteční hodnoty konstanty a náhodnou složku budou obsahovat tedy pouze přírůstky dX, poté pro kovarianční matici parametrů bude platit (s uvážením )

(37) .

(38) .

Budou-li parametry X vzájemně nezávislé, poté bude kovarianční matice diagonální s variancemi parametrů na diagonále. Zákon přenášení variancí (38) poté můžeme psát jako

(39) .

Odhad nejistoty parametrů vypočtených pomocí MNČ

(40) ,

kde je kovarianční matice měření. Podíváme-li se na vztah (18), ten je nezávislý na násobení váhové matice konstantou. Můžeme tak obecněji psát

(41) ,

kde je tzv. apriorní jednotková směrodatná variance (variance jednotkové váhy). Lze ukázat, že pro ní platí vztah

(42) ,

Dosazením (41) do (18) bychom dostali stejný výsledek. Konstantu tak můžeme volit, aniž bychom teoreticky ovlivnili výsledek vyrovnání.

Není-li předem apriorně známa, lze ji odhadnout z vyrovnaných měření. Lze ukázat, že tzv. nestranný odhad aposteriorní jednotkové variance je dán vztahem

(43) .

Zabývejme se nyní určením odhadu kovarianční matice vyrovnaných parametrů dle vztahu (18). Vztah (18) můžeme přeformulovat jako

(44) ,

kde a je někdy nazývána tzv. maticí kofaktorů. Užitím zákona přenášení variancí (32) poté dostáváme

(45)

(46)

(47) .

Příklad

Shrnutí odhadu parametrů metodou nejmenších čtverců a aposteriorního odhadu jejich nejistot

1. Problém vhodně formulujeme tak, aby šel maticově zapsat jako: .
2. Sestavíme nebo využijeme znalost kovarianční matice měření l k sestavení matice vah P: , nebo , kde použijeme buď apriorní jednotkovou varianci nebo ji nejprve volíme rovnu jedné, vypočteme její aposteriorní odhad a sestavení vah opakujeme. V případě měření stejné váhy volíme P jako jednotkovou matici.
3. Vypočteme vyrovnané parametry dle vztahu: .
4. Vypočteme odhad aposteriorní jednotkové variance: , kde N značí počet měření l a M počet parametrů q.
5. Vypočteme aposteriorní odhad kovarianční matice vyrovnaných parametrů: .
6. Směrodatné odchylky jednotlivých parametrů jsou dány odmocninou z diagonálních prvků kovarianční matice.
7. Případně vypočteme aposteriorní odhad kovarianční matice měření: .

Příklady

Aproximace polynomem

clear variables;

% Simulace dat polynomu

x = ( 0:5:50 )';

% Nominalni koeficienty

a_0 = 5; a_1 = 0.2; a_2 = 0.005; a_3 = -0.0001;

% Data y se zavedenim sumu

p = [ a_3, a_2, a_1, a_0 ];

y = polyval( p, x ) + 0.2*randn(size(x));

% Sestaveni matice A pro MNC, A_i = [ x^M, x^(M-1), ..., x^2, x, 1 ]; l je

% totozna s y

N = length(x);

A = ones( N, length(p) );

if length(p) > 1

for iN = 1:length(p)-1

A(:,iN) = x.^( length(p) - iN );

end

end

% Vypocet MNC

q = A\y;

% Residua

E = A*q - y;

% Nejistota vypoctenych parametru

sig_0_2 = E'*E/( length(y) - length(q) );

SIG_q = sig_0_2*inv( A'*A );

sig_q = sqrt( diag( SIG_q ) );

% Aposteriorni odhad nejistoty mereni

SIG_l = A*SIG_q*A';

sig_l = sqrt( diag( SIG_l ) );

% Tisk grafu s jemnejsim delenim aproximovanych dat

x_plot = linspace( min(x), max(x), 100 );

y_plot = polyval( q, x_plot );

figure

errorbar( x, y, sig_l, 'x' ) % Merena data i s aposteriornim odhadem nejistot mereni

hold on

plot( x_plot, y_plot, 'r-' ) % Aproximovana data

hold off

xlabel( 'x' )

ylabel( 'y' )

title_text = cell( length(p), 1 ); % Tvorba titulku formou bunky (vice radku popisku)

for it = 1:length(p)

title_text{it} = [ 'a_', num2str(length(p)-it), ' = ', num2str( q(it), '%.5f' ), ' (', num2str( sig_q( it ), '%.5f' ), ')' ];

end

title( title_text )

grid on

legend( 'měřená data \pm \sigma_{aposteriori}', 'aproximace', 'location', 'southeast' ) % Zobrazeni legendy

Aproximace přímých měření různé přesnosti

Vzdálenost L je měřena různými způsoby s různými nejistotami . Metodou MNČ vypočteme odhad skutečné hodnoty vzdálenosti a její nejistoty. Uvážíme-li nejistoty, poté výpočet MNČ odpovídá váženému průměru. Nebudeme-li uvažovat nejistoty, výpočet je poté roven průměru aritmetickému. Pro přehled jen uveďme použité vztahy pro aritmetický a vážený průměr, které lze snadno získat aplikací MNČ.

Aritmetický průměr

,

.

Vážený průměr

,

.

Aproximace plochy sférou

clear variables;

% Simulace dat na sfere

x = 0:5:50; y = x; % rozsah os

R = 80; % polomer nominalni sfery

x0 = 25; y0 = 25; % stred nominalni sfery

c = -1/R; % vrcholova krivost

[ X, Y ] = meshgrid( x, y ); % grid souradnic

RHO = sqrt( ( X - x0 ).^2 + (Y - y0).^2 ); % decentrace

% simulovana plocha i s sumem

Z = 25 + c*RHO.^2./( 1 + sqrt( 1 - c^2*RHO.^2 ) ) + 0.1*randn( size( X ) );

figure

mesh( X, Y, Z )

axis equal

xlabel( 'x' ); ylabel( 'y' ); zlabel( 'z' ); title( 'Měřená data' );

% Preskupeni matic do sloupcovych vektoru

[ M, N ] = size( X );

X = reshape( X, M*N, 1 );

Y = reshape( Y, M*N, 1 );

Z = reshape( Z, M*N, 1 );

% Sestaveni rovnic pro MNC

NN = length( X );

A = zeros( NN*(NN-1)/2, 3 );

l = zeros( NN*(NN-1)/2, 1 );

% Ak = [];

% lk = [];

idx1 = 0; idx2 = 0;

for iN = 1:NN-1

idx1 = idx2 + 1;

idx2 = idx1 + NN - iN - 1;

A( idx1:idx2, : ) = [ X(iN) - X(iN+1:end), Y(iN) - Y(iN+1:end), Z(iN) - Z(iN+1:end) ];

l( idx1:idx2 ) = 0.5*( ( X(iN)^2 + Y(iN)^2 + Z(iN)^2 ) - ( X(iN+1:end).^2 + Y(iN+1:end).^2 + Z(iN+1:end).^2 ) );

% Ak = [ Ak; X(iN) - X(iN+1:end), Y(iN) - Y(iN+1:end), Z(iN) - Z(iN+1:end) ];

% lk = [ lk; 0.5*( ( X(iN)^2 + Y(iN)^2 + Z(iN)^2 ) - ( X(iN+1:end).^2 + Y(iN+1:end).^2 + Z(iN+1:end).^2 ) ) ];

end

% Poloha stredu vypoctena pomoci MNC

rc = A\l;

% Vektor residui

E = A*rc - l;

% Nejistota polohy stredu

sig_0_2 = E'*E/( length( l ) - length( rc ) );

SIG_rc = sig_0_2*inv( A'*A );

sig_rc = sqrt( diag( SIG_rc ) );

disp( [ 'x_c = ', num2str( rc(1), '%.4f'), ' (', num2str( sig_rc(1), '%.4f') ')' ] );

x_c = 25.0313 (0.0056)

disp( [ 'y_c = ', num2str( rc(2), '%.4f'), ' (', num2str( sig_rc(2), '%.4f') ')' ] );

y_c = 25.1191 (0.0056)

disp( [ 'z_c = ', num2str( rc(3), '%.4f'), ' (', num2str( sig_rc(3), '%.4f') ')' ] );

z_c = -55.0052 (0.0433)

% Odhad polomeru sfery

Ri = sqrt( ( X - rc(1) ).^2 + ( Y - rc(2) ).^2 + ( Z - rc(3) ).^2 );

Rmnc = mean( Ri );

% Odchylky ve smeru polomeru

ER = Ri - Rmnc;

% Stredni kvadraticka odchylka polomeru

RMS = sqrt( ER'*ER/length( ER ));

% Histogram residui s hodnocenim RMS

figure

histogram( ER, 15 )

xlabel('\DeltaR')

ylabel('četnost')

title( { 'Histogram residuí \DeltaR' , [ 'RMS_R = ', num2str( RMS, '%.4f' ) ] } )