plik

6. Dopasowywanie modeli do danych

W procesie dopasowywania modelu matematycznego y

= g

(u,θ) do dodanych

eksperymentalnych wyodrębnić można dwa podstawowe etapy:

• identyfikację strukturalną, czyli dobór g

– odpowiedniego do problemu typu równań

matematycznych i powiązań między nimi (może się ona odbywać na drodze analizy
struktury fizycznej obiektu lub empirycznie);

• estymację parametrów, polegająca na doborze takich wartości parametrów modelu θ , aby

jak najdokładniej opisywał on posiadane dane (poza przypadkiem stosowania praw lub
teorii fizycznych, w których wartości parametrów są określone, ma ona charakter
empiryczny).

6.1. Estymacja parametrów jako zadanie odwrotne

Estymacja  parametrów  modelu  jest  zadaniem  odwrotnym  drugiego  rodzaju.  Na  podstawie
skutku,  jakim  są  dane  wyjściowe  y  (uzyskane  dla  pewnego  pobudzenia  u)  wyznaczana  jest
przyczyna takiej, a nie innej ich postaci, a są nią wartości parametrów modelu θ. Rozwiązanie
tego  zadania  wymaga  formalnego  odwrócenia  związku  przyczynowo-skutkowego
wyrażonego przez model g

, jego struktura mus być zatem znana.

Estymacja parametrów modeli fizykomatematycznych, gdzie odzwierciedlają one
właściwości fizyczne modelowanego obiektu, jest równoważna z pośrednim pomiarem
złożonym tych właściwości.

Można się zastanawiać, czy estymacja parametrów danego modelu (o znanej strukturze g

) jest jednakowo

efektywna dla dowolnego pobudzenia u, czy też istnieje pewne u optymalne, dające najdokładniejsze
oszacowanie θ. Na pytanie to odpowiada teoria planowania eksperymentów (rozdz. 8).

6.2. Estymacja parametrów modeli liniowych

Model może być liniowy np. względem wejść, parametrów, zmiennych stanu itd. W rozdziale
tym pod pojęciem „model liniowy” rozumiane będą modele liniowe względem parametrów
(LP), gdyż właśnie ta cecha determinuje postać algorytmów estymacji.

Przypomnijmy, że model liniowy musi spełniać zasadę superpozycji:





















gdzie a i b są skalarami, a całkowity wektor parametrów to θ





T
2



. Model taki,

uwzględniając przekształcenie pobudzenia u przez funkcję (liniową lub nieliniową) f

, można

ogólnie zapisać jako:

 



Techniki eksperymentu

lub posługując się macierzą eksperymentu (tożsamą z macierzą wrażliwości):

Xθ



gdzie dla N próbek i p parametrów:

 





























Zwykle liczba obserwacji (danych w spróbkowanych sygnałach u i y) jest większa od liczby
parametrów θ, tj. N > p.

Przykładem może być model w postaci szeregu potęgowego 3-go rzędu:

 











Wtedy macierz eksperymentu przyjmie postać:

























Ostatecznie zarejestrowane dane wyjściowe y to (w uproszczeniu) wyjście modelu y

zakłócone addytywnym szumem v:





6.2.1. Podstawowe założenia

Podczas estymacji parametrów modelu w oparciu o dane eksperymentalne najczęściej zakłada
się, że można je traktować jako realizację zmiennej losowej y opisanej funkcją gęstości
prawdopodobieństwa p(y,θ), co wypływa ze stochastycznego charakteru wektora parametrów

θ :

 



lub szumu v:









4. Wsadowe metody estymacji parametrów modeli liniowych

Prawdziwy (hipotetyczny) model związku przyczynowo-skutkowego o charakterze
deterministycznym oznaczony jest tu jako g

, a θ

to prawdziwe wartości parametrów tego

modelu.

Wektor v traktowany jest jako realizacja zmiennej losowej v. Jeżeli szum jest skorelowany, to
macierz jego kowariancji wynosi

 

 

νν





cov

W szczególności, gdy szum jest nieskorelowany, macierz kowariancji ma postać:









diag







lub dodatkowo dla szumu stacjonarnego:





Często zakłada się też, że addytywne zakłócenia losowe mają charakter szumu białego
(gaussowskiego):



, gdzie e to realizacja szumu białego o wariancji



Koncepcja wyjaśniająca losowy charakter danych wynikający z ich zakłócenia addytywnym
szumem pozwala następująco scharakteryzować właściwości statystyczne y:

 





 







cov

Natomiast w przypadku losowego charakteru wektora parametrów, θ

należy interpretować

jako

 

6.2.2. Metody wsadowe

Metody wsadowe pozwalają na obliczenie estymatorów parametrów modelu matematycznego
na podstawie wcześniej zarejestrowanych danych pomiarowych (off-line). W przypadku
modeli liniowych realizujące je algorytmy są jednokrokowe i stosunkowo nieskomplikowane.
Wadą ich jednak jest to, że estymacja może być przeprowadzona dopiero po zakończeniu
eksperymentu, a nie w czasie jego trwania.

6.2.2.1. Metody bayesowskie

Metody  bayesowskie  wykorzystują  prawdopodobieństwa  warunkowe.  Do  ich  grona  należy
metoda  maksymalnego  prawdopodobieństwa  a  posteriori  (MAP).  Dane  eksperymentalne  y
traktowane są w niej jako znane, natomiast  θ jako wektor losowy. Estymator

MAP

θˆ

uzyskuje

się maksymalizując prawdopodobieństwo tego, że pochodzi on z systemu o zarejestrowanym
(ustalonym) y:





max

arg

MAP



Techniki eksperymentu

gdzie funkcja gęstości prawdopodobieństwa a posteriori

 

θ |

na podstawie wzoru Bayesa

wynosi:







  

 



Funkcja gęstości prawdopodobieństwa

 

y |

znajdująca się w liczniuku nazywana jest

funkcją

wiarygodności

L(y,θ).

Ponieważ

jest

ustalone,

maksymalizacja

prawdopodobieństwa a posteriori sprowadza się do maksymalizacji licznika:

   

MAP

max

arg



W celu wyznaczenia wartości estymatora

MAP

θˆ

konieczna jest zatem aprioryczna znajomość

funkcji gęstości prawdopodobieństwa wektora parametrów f(



), co jednak rzadko ma miejsce.

W sprzyjających okolicznościach dostępne są typy poszczególnych rozkładów (rys.) oraz
estymaty ich dwóch pierwszych momentów, tj. wartości oczekiwanej i wariancji [Nied’09].

Rys. Przykładowe rozkłady gęstości prawdopodobieństwa parametru: normalny (lewy panel) i jednostajny

(prawy panel)

6.2.2.2. Metoda największej wiarygodności

W metodzie największej wiarygodności (ML) dane eksperymentalne y traktowane są jako
zmienna losowa, a



jako nieznany, lecz ustalony wektor parametrów. Estymację można

interpretować  jako  wybór  takiego  wektora  parametrów,  dla  którego  zarejestrowane  dane  są
najbardziej prawdopodobne (wg funkcji wiarygodności L).

Przy maksymalizacji prawdopodobieństwa wystąpienia zaobserwowanych wartości y (tzn. ich
wiarygodności)  postępuje  się  analogicznie  jak  przy  dopasowaniu  odpowiedzi  modelu  do
pomiarów tyle, że tym razem model odpowiedzi ma charakter statystyczny, podany w postaci
funkcji gęstości prawdopodobieństwa wystąpienia ich określonych wartości.

Funkcja gęstości prawdopodobieństwa ciągu próbek

y ,...,

sygnału odpowiedzi na znane

pobudzenie ma postać wielowymiarową, a wartości próbek zależą od poszukiwanych
parametrów modelu θ. Każda z próbek jest realizacją pewnej zmiennej losowej y

, a funkcja

gęstości opisuje cały wektor zmiennych





y ,...,



w sposób łączny.

4. Wsadowe metody estymacji parametrów modeli liniowych

Estymator

θˆ

maksymalizuje funkcję wiarygodności L(y,θ) ≡ f(y|θ):

 

max

arg



Maksymalizacja może być przeprowadzona tylko przy znajomości postaci funkcji
wiarygodności. Zatem, aby obliczyć

θˆ

, należy znać lub założyć na podstawie posiadanej

wiedzy jej jawną postać zależną od θ.

W  ogólności  poszukiwanie  ekstremum  nieliniowej  funkcji  L  sprowadza  się  do  problemu
optymalizacyjnego  (patrz  rozdz.  3.4).  W  przypadku  pewnych  rozkładów  gęstości
prawdopodobieństwa  zadanie  upraszcza  się  jednak  do  estymacji  liniowej.  Takim
przypadkiem,  najczęściej  też  stosowanym  w  praktyce,  jest  opis  addytywnych  zakłóceń
pomiaru  poprzez  rozkład  normalny,  a  w  szczególności  przez  ich  macierz  kowariancji  R.
Wtedy  wielowymiarowy  rozkład  łączny  ma  funkcję  gęstości  prawdopodobieństwa  (funkcję
wiarygodności) postaci:



 

  

 





 





exp

det







gdzie wektor rezyduów (reszt):

ε(θ) = y – y

(θ).

Korzystając z faktu, iż logarytm jest funkcją monotoniczną, poszukiwanie maksimum funkcji
wiarygodności sprowadza się do poszukiwaniu maksimum jej logarytmu naturalnego postaci:

 





 

det





co upraszcza obliczenia. Ponieważ pierwsze dwa składniki maksymalizowanego wyrażenia
nie zależą od estymowanych parametrów, to ostatecznie maksymalizacja funkcji
wiarygodności prowadzi do minimalizacji (znak minus przy ostatnim czynniku) funkcjonału

 





 





 





min

arg

max

arg







Różniczkując macierzowo funkcję kryterialną V

(θ) względem θ

 





 

























Xθ





















i przyrównując pochodną do zera:

















Techniki eksperymentu

Znajduje się minimum V

dla estymatora postaci:









o macierzy kowariancji S

 





cov





Warto zauważyć, że estymator ML jest szczególnym przypadkiem estymatora MAP wtedy,
gdy rozkłady gęstości prawdopodobieństwa parametrów modelu są jednostajne:

 

MAP

cons



gdyż wtedy:

   

 

MAP

max

arg

max

arg



6.2.2.3. Metoda najmniejszych kwadratów

W metodzie najmniejszych kwadratów (LS), wprowadzonej przez Legendre’a i Gaussa,
zakłada się addytywny charakter szumu będącego realizacją nieskorelowanych zmiennych
losowych o rozkładzie normalnym, zerowej wartości oczekiwanej i jednakowej wariancji (co
w wielu wypadkach jest uzasadnione dzięki centralnemu twierdzeniu granicznemu). Podejście
to wymaga zatem niewielkiej wiedzy apriorycznej.

Estymator

θˆ uzyskuje się minimalizując odległość modelu y

(θ) od danych

eksperymentalnych y w przestrzeni l

(odległość euklidesowa):

 

min

arg



gdzie funkcjonał V

ma postać:

 











 



Xθ





















Zauważmy, że podejście to stosowane jest przy założeniu, że addytywne błędy pomiarowe e
występują tylko w zmiennej zależnej (zarejestrowany sygnał wyjściowy y), natomiast
zmienna niezależna (pobudzenie u) znana jest bezbłędnie (rys).

4. Wsadowe metody estymacji parametrów modeli liniowych

Rys. Ilustracja błędów zmiennej zależnej y w metodzie najmniejszych kwadratów

W celu znalezienia wektora parametrów dającego minimum funkcjonału należy
zróżniczkować V

względem θ:





Xθ

)

(













i uzyskane wyrażenie na pochodną przyrównać do zera:







Macierz X

X o wymiarze p



p (gdzie p to liczba parametrów) jest kwadratowa i jeżeli jest

nieosobliwa, to można ją odwrócić uzyskując wyrażenie na estymator:

 





Warunkiem na to, że wyznaczone ekstremum to minimum, jest dodatnio określony hessian,
czyli macierz drugich pochodnych:

 





Xθ













Macierz X

X jest dodatnio określona, zatem znalezione ekstremum to minimum.

Macierz kowariancji estymatora LS dana jest przez:

 

cov







gdzie



, dla przypomnienia, jest wariancją addytywnego szumu białego. Ponieważ metodę

LS stosuje się w sytuacjach, w których nieznane są właściwości statystyczne zakłóceń, nie
można wyznaczyć odchyleń standardowych estymatorów bezpośrednio z powyższego wzoru
(wartość



nie jest znana). Wariancję szumu można jednak oszacować na podstawie

Techniki eksperymentu

rezyduów, zakładając bezbłędność modelu. Wtedy błąd średniokwadratowy V

MSE

dopasowania modelu do danych wynosi [Sjöberg’95]:

 





MSE









Pamiętając, że jednocześnie

 





 









otrzymuje się znane wyrażenie pozwalające oszacować wariancję szumu na podstawie
wyznaczonych rezyduów:







i ostatecznie [Hocking’76]:

   

 

















Estymator LS jest szczególnym przypadkiem estymatora ML (a tym samym MAP) w sytuacji,
gdy addytywne zakłócenia mają charakter stacjonarny (tj.









) i są

realizacją niezależnych zmiennych losowych (tzn. macierz kowariancji jest diagonalna):





gdyż wtedy:









 







Przykład: Estymacja parametrów układu oddechowego podczas sztucznej wentylacji

W uproszczeniu układ oddechowy można zamodelować jako szeregowe połączenie dwóch parametrów: oporu
R

i podatności C

(rezystor i kondensator w elektrycznym modelu zastępczym – rys.).

Rys. Elektryczny model zastępczy procesu sztucznej wentylacji płuc

Z II prawa Kirchhoffa otrzymuje się związek (model MISO) między wyjściem (ciśnienie P

mierzone w

respiratorze) i wejściem modelu (przepływ Q mierzony w respiratorze oraz obliczana na jego podstawie objętość
powietrza wprowadzonego do płuc V = ∫Qdt ):

 





4. Wsadowe metody estymacji parametrów modeli liniowych

który w zapisie macierzowym przyjmuje postać (przez P

oznaczając ciśnienie na początku pomiarów):





Wektor parametrów modelu to θ = [R

, E

, P

]

. Macierz eksperymentu wynosi:











zatem estymaty parametrów (w tym mechanicznych właściwości układu oddechowego) oblicza się następująco:













Przykładowy przepływ wentylacyjny Q oraz wynik dopasowania modelu do danych ciśnienia P

pokazano na

rys.

Rys. Przykładowe dopasowanie modelu sztucznej wentylacji płuc do danych

6.2.2.3. Metoda ważonych najmniejszych kwadratów

Często stosowanym, uogólnionym podejściem do estymacji, jest metoda ważonych
najmniejszych kwadratów (WLS), w której estymatory

WLS

θˆ

wyznacza się minimalizując

funkcjonał V

WLS

 

WLS

min

arg



w postaci ważonego błędu średniokwadratowego (odległości modelu od danych w l

 





 



 



Xθ









WLS

gdzie W jest różnie w poszczególnych metodach dobieraną, dodatnio określoną, symetryczną
macierzą wag (najczęściej diagonalną). Wyrażają one posiadaną wiedzę aprioryczną o
poszczególnych danych pomiarowych, np. ich dokładność czy szczególne znaczenie.

Wektor estymatorów i macierz kowariancji dane są w tym przypadku wzorami [VBos’88]:

0.5

1.5

2.5

3.5

-1

-0.5

0.5

Model 1

(

)

0.5

1.5

2.5

3.5

-0.1

0.1

0.2

0.3

0.4

Time (s)

(

)

Data
Model

Techniki eksperymentu









WLS

 









cov





WRWX

WLS

gdzie R, przypomnijmy, jest dodatnio określoną macierzą kowariancji zakłóceń addytywnych
o zerowej wartości oczekiwanej.

Nietrudno zauważyć, że zarówno ML jak i LS są szczególnymi przypadkami metody
ważonych najmniejszych kwadratów. Można to sprawdzić podstawiając W = R

–1

w przypadku

metody największej wiarygodności (wagi są wtedy odwrotnie proporcjonalne do niepewności
poszczególnych  danych)  lub  W = I  w  przypadku  metody  najmniejszych  kwadratów
(wszystkie wagi są jednakowe i wynoszą 1).

Zależność estymatora od wyjścia modelu i szumu pomiarowego
Na  przykładzie  metody  WLS  warto  przyjrzeć  się,  w  jaki  sposób  estymator  parametrów
budowany  jest  przez  wyjście  modelu  y

(tożsame z niezakłóconym wyjściem obiektu w

przypadku bezbłędnego modelu), a w jakim stopniu przez konkretną realizację szumu
pomiarowego e (pamiętając, że w przypadku szumu addytywnego obowiązuje relacja
y = y

+ e):

































WXθ



















WLS

gdzie θ

oznacza wektor parametrów o prawdziwych wartościach (tj. takich, które miał model

generujący wyjście y

= Xθ

). Widać stąd, że wpływ obu czynników (y

i e) jest analogiczny

– są one tak samo przetworzone przez operator







. W efekcie estymator

WLS

θˆ

ma charakter zmiennej losowej (także w innych metodach), a obliczane jego wartości tym
bardziej są zbliżone do wartości oczekiwanej, im mniejsza jest energia zakłóceń losowych
zawartych w danych empirycznych (im większy stosunek sygnału do szumu – SNR).

6.2.3. Estymacja z ograniczeniami

6.2.3.1. Estymacja z ograniczeniami liniowymi

W  rzeczywistych  sytuacjach  parametry  odpowiadające  fizycznym  właściwościom  obiektu
mogą  przyjmować  tylko  pewne  wartości,  zazwyczaj  należące  do  znanego  przedziału  [a, b],
lub  np.  nie  mogą  być  ujemne.  W  procesie  wyznaczania  ich  wartości  należy  wykorzystać  tę
informację aprioryczną, a stosowane metody noszą nazwę estymacji z ograniczeniami.

4. Wsadowe metody estymacji parametrów modeli liniowych

Zdanie estymacji z ograniczeniami jest złożone i zaproponowano wiele metod jego
rozwiązania. Do najpopularniejszych należą metody liniowe (opierające się o LS), w których
ograniczenia nałożone na θ przyjmują postać układu równań liniowych:

Gθ



gdzie G

(q×p)

(q < p) jest znaną macierzą więzów, a g jest znanym wektorem współczynników

liniowych. Należy zatem rozwiązać typowe zadanie LS ale tak, aby estymatory parametrów
spełniały powyższe równanie:



 

















min

arg

Xθ

W  tym  celu  stosuje  się  jedną  z  dwóch  metod:  metodę  mnożników  Lagrange’a  (typowo  gdy
q < p/2) lub metodę elementów (q > p/2).

Metoda elementów
Metoda  elementów  wykorzystuje  fakt,  że  skoro  p  parametrów  jest  skojarzonych  ze  sobą  za
pośrednictwem q równań liniowych, to niezależnych jest tylko q z nich, a inne można wyrazić
jako  kombinacje  liniowe  pozostałych.  Ostatecznie  wystarczy  wyestymować  wartości  p – q
parametrów,  a  q  następnych  obliczyć  na  ich  podstawie.  Prowadzi  to  do  redukcji  liczby
estymowanych parametrów i stąd metoda ta zwana jest również metodą redukcji zmiennych.

Pierwszym krokiem jest przedstawienie macierzy G jako macierzy blokowej, podzielonej na
macierz  prostokątną  G

((p–q)×q)

oraz macierz kwadratową G

(q×q)

: G [G

], a wektora

parametrów odpowiednio jako [θ

]

. Wtedy układ równań więzów przyjmuje postać:



















Stąd:





a rozwiązanie względem θ











Pozwala to na wyeliminowanie q parametrów z równania modelowego, przy czym macierz
eksperymentu również dzieli się na odpowiednie bloki:























Podstawiając teraz wyprowadzoną zależność na θ

, otrzymuje się następujące równanie z

niewiadomym wektorem p – q parametrów:













Techniki eksperymentu

Można je uporządkować do podstawowej postaci liniowej (przenosząc na lewo znane
wartości) i pamiętając, że posiadane dane y = y

+ e, zapisać:













 



 









a następnie znaleźć estymator θ

metodą najmniejszych kwadratów. Ostatecznie uzyskuje się:























6.2.4. Estymacja parametrów modeli fourierowskich

Analiza  właściwości  danych  powiązanych  z  ich  widmem  nie  musi  być  przeprowadzana  w
dziedzinie  częstotliwości.  W  wielu  sytuacjach  znane  są  częstotliwości  składowe
analizowanych  sygnałów,  a  poszukiwane  jedynie  ich  amplitudy  i  przesunięcia  fazowe.
Przydatne  okazuje  się  wtedy  wykorzystanie  modelu  matematycznego  badanego  sygnału
analogicznego do szeregu Fouriera (częstotliwości nie muszą być harmoniczne):

 

























sin

cos

sin







gdzie:







sin







Ponieważ w tym przypadku nieznanymi parametrami są α

, a i b, dopasowywanie tego

modelu do danych sprowadza się do estymacji modelu LP w dziedzinie czasu (rys.), a
macierz eksperymentu przyjmuje postać:































































sin

cos

sin

cos

sin

cos

sin

cos

sin

cos

sin

cos









