Rachunek prawdopodobieństwa i statystyka/Wykład 12: Metoda największej wiarygodności

Metoda największej wiarygodności

Omówimy metodę estymacji punktowej, zwaną metodą największej wiarygodności oraz pokażemy konkretne problemy, w których może być ona stosowana. Zwrócimy także uwagę na trudności, jakie można często spotkać stosując tę metodę. Najpierw jednak przypomnimy pewne fakty dotyczące optymalizacji funkcji.

Wartość największa funkcji

Poznamy teraz chyba najpopularniejszą metodę estymacji punktowej - metodę największej wiarygodności. Jednak aby ją poprawnie stosować, musimy przypomnieć sobie pewne wiadomości z analizy matematycznej (patrz wykład z Analizy matematycznej).

Przypuśćmy, że mamy daną funkcję $f\colon K\longrightarrow {\mathbb {R} }$, gdzie $K\subset {\mathbb {R} }$ jest ustalonym zbiorem. Mówimy, że funkcja ta przyjmuje wartość największą w punkcie ${\hat {x}}\in K$, jeżeli:

$f(x)\leq f({\hat {x}})\;\;$ dla każdego $\;x\in K$

Oczywiście, nie dla wszystkich funkcji daje się określić wartość największą, jednak przy pewnych dodatkowych założeniach można stwierdzić, że wartość taka istnieje. Mówi o tym poniższe twierdzenie, które przytaczamy bez dowodu.

Twierdzenie 12.1

Załóżmy, że funkcja $f$ jest ciągła na zbiorze $A$ oraz że zachodzi jeden z następujących warunków:

1. $A=[a,b]$ jest przedziałem domkniętym i ograniczonym,

2. $A$ jest dowolnym przedziałem (ograniczonym lub nieograniczonym) oraz istnieją granice funkcji $f$ na końcach tego przedziału i są one skończone lub równe $-\infty$.

Wtedy funkcja $f$ przyjmuje wartość największą w pewnym punkcie przedziału $A$.

Z praktycznego punktu widzenia, zainteresowani jesteśmy wyznaczeniem punktu ${\hat {x}}$, w którym dana funkcja przyjmuje wartość największą. Bardzo ważnym narzędziem okazuje się być tutaj pochodna (patrz wykład z Analizy matematycznej) - oto, bardzo pomocne w tym przypadku, klasyczne twierdzenie z analizy matematycznej:

Twierdzenie 12.2

Jeżeli funkcja $f\colon (a,b)\longrightarrow {\mathbb {R} }$ jest różniczkowalna i przyjmuje wartość największą w punkcie ${\hat {x}}\in (a,b)$, to $f'({\hat {x}})=0$.

Podkreślamy, że w obu powyższych twierdzeniach wszystkie założenia są istotne. Jeżeli w konkretnej sytuacji potrafimy stwierdzić, że są one spełnione, to nasz problem sprowadza się do obliczenia pochodnej i rozwiązania równania:

$f'(x)=0$

Wówczas funkcja $f$ może osiągać wartość największą jedynie w punktach będących rozwiązaniami powyższego równania lub końcami przedziału określoności, o ile należą one do tego przedziału. Bardzo często zdarza się, że nasze równanie ma dokładnie jeden pierwiastek oraz że łatwo sprawdzić, iż wartość największa nie może być przyjęta na końcach przedziału określoności - w tym przypadku to właśnie owo rozwiązanie jest jedynym punktem, w którym funkcja przyjmuje wartość największą.

W niektórych przypadkach funkcja $f$ jest na tyle skomplikowana, że nie potrafimy stwierdzić, czy zachodzą założenia twierdzenia 12.1 i twierdzenia 12.2. Praktyczną metodą jest wtedy narysowanie wykresu (na przykład za pomocą komputera) i na zauważenie na jego podstawie, że taka wartość rzeczywiście istnieje. Innym problemem może być brak różniczkowalności lub skomplikowana postać pochodnej $f'(x)$, uniemożliwiająca analityczne rozwiązanie powyższego równania - należy wtedy zastosować odpowiednią metodę numeryczną.

Podkreślamy, iż metoda największej wiarygodności, którą za chwilę przedstawimy, jest zaimplementowana w większości komputerowych programów matematycznych i statystycznych. Na przykład, program Maple (w wersji 10) udostępnia ją w pakiecie:

Statistics[MaximumLikelihoodEstimate],

zaś w programie Excel istnieje dodatek Solver, który można, między innymi, zastosować do optymalizacji funkcji.

Estymacja metodą największej wiarygodności

Omówimy tutaj jedną z najczęściej stosowanych metod estymacji punktowej - metodę największejwiarygodności. Zaczniemy od (fikcyjnego) przykładu.

Przykład 12.3

Spośród studentów informatyki pewnego elitarnego wydziału wybrano losowo i niezależnie od siebie 50 osób, a następnie każdą z nich spytano, czy kiedykolwiek w trakcie studiów otrzymała ocenę niedostateczną. Okazało się, iż 14 osób odpowiedziało "TAK", natomiast pozostałe odpowiedziały "NIE". Pytamy teraz: jaki procent studentów informatyki otrzymał w trakcie swoich studiów ocenę niedostateczną.

Mamy tutaj zaobserwowaną próbkę prostą $x_{1},\dots ,x_{n}$, $n=50$, z rozkładu dwupunktowego $(0,1,p)$: $0$ interpretujemy jako "NIE", zaś $1$ - jako "TAK". Naszym zadaniem jest wskazanie parametru $p$. Oczywiście, nie potrafimy tego zrobić dokładnie na podstawie samej tylko próbki, natomiast możemy możliwie najlepiej przybliżyć jego nieznaną wartość w następujący sposób: obliczamy prawdopodobieństwo wylosowania naszej próbki w zależności od $p$, a następnie uznajemy, że najlepszym przybliżeniem nieznanego parametru będzie taka wartość $p$, dla której obliczone właśnie prawdopodobieństwo jest największe.

Przystąpmy zatem do realizacji opisanej powyżej procedury. Korzystając z niezależności zmiennych losowych $X_{1},\dots ,X_{n}$ otrzymujemy:

$P(X_{1}=x_{1},\dots ,X_{n}=x_{n})=P(X_{1}=x_{1})\cdot \dots \cdot P(X_{n}=x_{n})$

Zauważmy, że:

$P(X_{i}=x_{i})=\left\{{\begin{array}{rl}p,&{\text{gdy}}x_{i}=1\\1-p,&{\text{gdy}}x_{i}=0\end{array}}\right.$

Z treści zadania wiemy, że $x_{i}=1$ dla dokładnie 14 wartości $i$ . Tak więc:

$P(X_{1}=x_{1},\dots ,X_{n}=x_{n})=p^{\sum x_{i}}(1-p)^{n-\sum x_{i}}=p^{14}(1-p)^{36}$

Pozostaje nam wyznaczyć największą wartość funkcji $l\colon [0,1]\longrightarrow {\mathbb {R} }$, zadanej wzorem:

$l(p)=p^{14}(1-p)^{36}$

oraz zwanej funkcją (największej) wiarygodności. Łatwo stwierdzić, że funkcja ta ma wartość największą, gdyż jest ciągła na przedziale domkniętym i ograniczonym $[0,1]$. Co więcej, wartość ta musi być przyjęta w jakimś punkcie (lub punktach) ${\hat {p}}\in (0,1)$, gdyż dla $p=0$ oraz dla Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle p= 1} wartości funkcji Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle l} są równe Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle 0} , i właśnie to Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle \hat{p}} przybliża nieznaną wartość parametru Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle p} .

W celu wyznaczenia Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle \hat{p}} wykorzystamy powszechnie używaną metodę upraszczającą obliczenia - rozważymy mianowicie funkcję:

Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle L(p) = \ln l(p) }

która przyjmuje wartość największą dokładnie w tych samych punktach, co funkcja Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle l} . Tak więc:

Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle L(p) = 14 \ln p + 36 \ln (1-p) }

Obliczamy:

Parser nie mógł rozpoznać (SVG (MathML może zostać włączone przez wtyczkę w przeglądarce): Nieprawidłowa odpowiedź („Math extension cannot connect to Restbase.”) z serwera „https://wazniak.mimuw.edu.pl/api/rest_v1/”:): {\displaystyle L'(p) = \frac{14}{p} - \frac{36}{1-p} }

a następnie rozwiązujemy równanie $L'(p)=0$, czyli:

${\frac {14}{p}}-{\frac {36}{1-p}}=0$

otrzymując następujące rozwiązanie:

${\hat {p}}={\frac {14}{50}}=0.28$

Otrzymany w ten sposób estymator nazywa się estymatorem największej wiarygodności parametru $p$.

Metoda największej wiarygodności polega więc na skonstruowaniu funkcji wiarygodności odpowiadającej zaobserwowanemu zdarzeniu, zależnej od szukanych (estymowanych) parametrów, a następnie na znalezieniu takich wartości tych parametrów, dla których funkcja ta osiąga największą wartość. Podkreślamy jednak, że wartość funkcji największej wiarygodności nie musi być dokładnie równa prawdopodobieństwu zaobserwowanego zdarzenia - wystarczy, że będzie do niego proporcjonalna (patrz przykład 12.4).

Przykład 12.4

Chcąc zbadać wadliwość nowej serii komputerów, przeprowadzono następujące badanie: przez 20 dni uruchamiano codziennie 10 nowych komputerów i każdy z nich poddawano wszechstronnemu testowi. Otrzymano następujące wyniki: w ciągu 14 dni wszystkie komputery działały bez zarzutu, w ciągu 4 dni miała miejsce awaria jednego z komputerów, natomiast w ciągu 2 dni zaobserwowano awarie 2 komputerów. Jaka jest wadliwość losowo wybranego

komputera, rozumiana jako prawdopodobieństwo awarii w czasie jednego dnia pracy?

Oznaczmy szukaną wadliwość komputera przez $p$ i policzmy prawdopodobieństwo zaobserwowanego zdarzenia w zależności od $p$. W tym celu zauważmy najpierw, że prawdopodobieństwo zajścia dokładnie $k$ awarii w ciągu jednego dnia wynosi:

$a_{k}={\dbinom {10}{k}}p^{k}(1-p)^{10-k}$

Ponieważ awarie zachodzą niezależnie od siebie, więc prawdopodobieństwo opisanego powyżej zdarzenia wyraża się wzorem:

$a_{0}^{14}a_{1}^{4}a_{2}^{2}=\left((1-p)^{10}\right)^{14}\left({\dbinom {10}{1}}p(1-p)^{9}\right)^{4}\left({\dbinom {10}{2}}p^{2}(1-p)^{8}\right)^{2}$

$=20250000(1-p)^{192}p^{8}$

Jako funkcję wiarygodności warto więc przyjąć:

$l(p)=(1-p)^{192}p^{8}$

Naszym zadaniem jest znalezienie takiego punktu ${\hat {p}}$, w którym funkcja $l$ osiąga wartość największą na przedziale $[0,1]$. Zauważmy, że taka wartość ${\hat {p}}$ istnieje i jest liczbą z przedziału $(0,1)$. Aby ją wyliczyć postępujemy dokładnie tak samo, jak poprzednio - definiujemy:

$L(p)=\ln l(p)=192\ln(1-p)+8\ln p$

obliczamy pochodną:

$L'(p)=-{\frac {192}{1-p}}+{\frac {8}{p}}$

a następnie rozwiązujemy równanie $L'(p)=0$, otrzymując:

${\hat {p}}=0.04$

Do tej pory rozważaliśmy jedynie przykładowe sytuacje, w których miała zastosowanie metoda największej wiarygodności. Zajmijmy się więc teraz przypadkiem ogólnym.

Jeżeli obserwujemy próbkę prostą $x_{1},\dots ,x_{n}$z rozkładu dyskretnego o parametrze $\theta$, to określamy funkcję wiarygodności jako:

$l(\theta )=cP_{\theta }(x_{1})\cdot \dots \cdot P_{\theta }(x_{n})$

gdzie $c$ jest stałą dodatnią, zaś estymatorem największej wiarygodności parametru $\theta$ nazywamy taką wartość ${\hat {\theta }}\in \Theta$, że dla każdego $\theta \in \Theta$ zachodzi warunek:

$l({\hat {\theta }})\leq l(\theta )$

W przypadku rozkładów ciągłych prawdopodobieństwo zaobserwowania pojedynczej próbki prostej $x_{1},\dots ,x_{n}$ jest równe $0$, jednak i w tym przypadku można stosować metodę największej wiarygodności - tym celu definiuje się funkcję wiarygodności:

$l(\theta )=cf_{\theta }(x_{1})\cdot \dots \cdot f_{\theta }(x_{n})$

gdzie $f_{\theta }$ jest gęstością rozkładu $P_{\theta }$, zaś $c>0$ jest stałą.

Przykład 12.5

Rozważmy próbkę prostą z rozkładu wykładniczego o parametrze $\lambda >0$. Znajdziemy

estymator największej wiarygodności dla tego parametru.

Pamiętamy z wykładu 8 że gęstością rozkładu wykładniczego jest funkcja:

$f(x)=\left\{{\begin{array}{rl}0&{\text{dla}}x<0\\\lambda e^{-\lambda x}&{\text{dla}}x\geq 0\end{array}}\right.$

Ponieważ próbka $x_{1},\dots ,x_{n}$ pochodzi z tego rozkładu, więc można założyć, że:

$x_{i}>0\;\;$ dla każdego $i$

Zatem funkcja wiarygodności ma w tym przypadku postać:

$l(\lambda )=\lambda e^{-\lambda x_{1}}\cdot \dots \cdot \lambda e^{-\lambda x_{n}}=\lambda ^{n}e^{-\lambda \sum _{i=1}^{n}x_{i}}=\lambda ^{n}e^{-\lambda n{\bar {x}}}$

(tutaj $n{\bar {x}}>0$ jest znaną liczbą). Jak widać, funkcja $l$ ma w punkcie $0$ wartość równą $0$, a także można łatwo stwierdzić, że:

$\lim _{\lambda \longrightarrow \infty }l(\lambda )=0$

Jest oczywiste, że $l$ jest funkcją ciągłą, przyjmującą wartości dodatnie dla wszystkich $\lambda >0$. Tak więc istnieje punkt ${\hat {\lambda }}>0$, w którym funkcja $l$ przyjmuje wartość największą. Aby go wyznaczyć, wygodnie jest rozważyć funkcję:

$L(\lambda )=\ln {l(\lambda )}=n\ln \lambda -\lambda n{\bar {x}}$

Teraz różniczkujemy:

$L'(\lambda )={\frac {n}{\lambda }}-n{\bar {x}}$

i widzimy, że pochodna $L'$ przyjmuje wartość zero w punkcie:

${\hat {\lambda }}={\frac {1}{\bar {x}}}$

który jest właśnie szukanym estymatorem parametru $\lambda$.

Estymatory największej wiarygodności - własności

Poznaliśmy ogólne zasady konstrukcji estymatorów metodą największej wiarygodności. Jednak uważny student zwrócił z pewnością uwagę, że w niektórych przypadkach można było z góry przewidzieć wynik. Czy musieliśmy więc używać wówczas metody? Oczywiście nie - metodę największej wiarygodności stosuje się, przede wszystkim, w sytuacjach, w których nie widać od razu rozsądnego estymatora (w trakcie ćwiczeń omówimy takie sytuacje). Jednakże warto zwrócić uwagę na to, że estymatory największej wiarygodności posiadają pewne uniwersalne własności, co sprawia, że są one na ogół "dobrymi" estymatorami. Poniżej przytaczamy niektóre z tych własności. Pamiętajmy jednak, iż, aby one zachodziły, należy przyjąć pewne dość techniczne założenia, które na ogół są spełnione.

(1) Estymator największej wiarygodności jest zgodny.

(2) Estymator największej wiarygodności jest asymptotycznie nieobciążony.

(3) W przypadku dużych próbek, estymator największej wiarygodności parametru $\theta$ ma w przybliżeniu rozkład $N(\theta ,{\frac {1}{\sqrt {I_{n}}}})$, gdzie $I_{n}$ jest tak zwaną informacją Fishera, którą można określić dla niemal każdego rozkładu dyskretnego lub ciągłego (nie robimy tego jednak tutaj).

(4) Jeżeli ${\hat {\theta }}$ jest estymatorem największej wiarygodności parametru $\theta$, zaś $g\colon \Theta \longrightarrow {\mathbb {R} }$ - funkcją ciągłą, to $g({\hat {\theta }})$ jest estymatorem największej wiarygodności parametru $g(\theta )$.

Metodę największej wiarygodności stosuje się także w sytuacji, gdy szukany parametr jest wektorem, na przykład $\theta =(m,\sigma )$ w rozkładzie $N(m,\sigma )$. Należy wówczas wyznaczyć wartość największą funkcji wielu zmiennych (patrz wykład z Analizy matematycznej 2), co jednak często okazuje się być zadaniem niezbyt łatwym.