Naszym celem jest przewidywanie przyszłych szkód na podstawie zgroma-dzonych przez ubezpieczyciela danych. W modelu bayesowskim, zadanie predykcji ma eleganckie rozwiązanie, które przedstawiliśmy w poprzednim podrozdziale. Niestety, użyteczność tego rozwiązania jest ograniczona. Obli-czenie najlepszego predyktora wymaga znajomości łącznego rozkładu praw-dopodobieństwa rozważanych zmiennych losowych. Dlatego warto zadanie postawić trochę mniej ambitnie i skupić uwagę na predykcji liniowej. Dzięki temu można rozważać modele znacznie prostsze i wymagające mniej szcze-gółowych założeń.
Predykcja liniowa
Rozważamy, podobnie jak w poprzednim podrozdziale, wektor obserwacji X = (X1, . . . , Xn) i nieobserwowalną zmienną losową Y . Poszukujemy teraz predyktora h(X) zmiennej Y , który jest liniową funkcją obserwacji:
h(X) = h(X1, . . . , Xn) = c0 +
n
X
i=1
ciXi.
W dalszym ciągu, mówiąc o funkcjach liniowych mamy na myśli funkcje nie-jednorodne, które mogą zawierać wyraz wolny. Tak jak poprzednio, za kry-terium jakości przyjmiemy błąd średniokwadratowy MSE. Od tego momentu zakładamy stale, że EY2 < ∞ i EXi2 < ∞.
4.4.1 DEFINICJA. Zmienna losowa ˆY = ˆh(X) jest najlepszym liniowym predyktorem zmiennej Y , jeśli ˆh : Rn→ R jest funkcją liniową i
E Y −ˆh(X)2
≤ E Y − h(X)2
dla każdej funkcji liniowej h : Rn → R. Będziemy symbolicznie pisali ˆY =
BLP(Y ).
Skrót BLP pochodzi od „Best Linear Predictor”.
4.4.2 TWIERDZENIE (Najlepszy liniowy predyktor). Zmienna losowa c0+Pn
i=1ciXi jest najlepszym liniowym predyktorem Y , jeśli współczynniki c0, c1, . . . , cn spełniają układ równań:
n
X
i=1
ciCov(Xi, Xk) = Cov(Xk, Y ), k = 1, . . . , n;
c0 = EY −
n
X
i=1
ciEXi.
Dowód. Ponieważ ograniczamy się do funkcji liniowych, MSE jest faktycznie
Czytelnik łatwo przekona się, że te iloczyny mieszane w rozwinięciu kwa-dratu, które nie zostały napisane, mają wartość oczekiwaną zero.
Tak więc, MSE =
Pierwszy składnik w ogóle nie zależy od współczynników ci. Z kolei ostatni składnik można zawsze zmniejszyć do zera przez odpowienie dobranie wyrazu wolnego c0, dla każdych c1, . . . , cn. Ostatnie z równań w tezie twierdzenia dostajemy, przyrównując do zera ten składnik. Pozostaje więc wyznaczyć c1, . . . , cn tak, żeby zminimalizować dwa środkowe składniki. Mamy
Przyrównując pochodne cząstkowe do zera, otrzymujemy nasz układ równań.
Układ równań w tym twierdzeniu zawsze ma rozwiązanie. To po prostu dla-tego, że nieujemna funkcja kwadratowa gdzieś musi przyjmować najmniejszą wartość. Niekiedy rozwiązań może być dużo. Najczęściej będziemy mieli do czynienia z „niezdegenerowanym przypadkiem”, kiedy nasz układ n + 1 równań z n + 1 niewiadomymi ma jedno rozwiązanie.
Ważne dla nas będzie spostrzeżenie, że możemy wyznaczyć najlepszy liniowy predyktor bez dokładnej znajomości rozkładów prawdopodobieństwa. Wy-starczy znać wartości oczekiwane i kowariancje wszystkich interesujących nas zmiennych.
Liniowa predykcja w modelu bayesowskim
Wróćmy do modelu bayesowskiego z poprzedniego rozdziału. Rozpatrujemy zmienne losowe Θ, X1, . . . , Xnspełniające Założenie 4.3.3. W tym modelu, Θ jest losową „zmienną strukturalną”, z samej swojej natury nieobserwowalną.
Interesuje nas predykcja µ(Θ) na podstawie obserwacji X1, . . . , Xn.
Wartości oczekiwane i struktura kowariancji. Przypomnijmy, że funk-cja µ jest określona następująco.
µ(θ) = Eθ(Xi) =
(R xfθ(x)dx dla zmiennej typu ciągłego;
P
xxfθ(x) dla zmiennej typu dyskretnego.
Analogicznie, niech
σ2(θ) = Varθ(Xi) =
(R(x − µ(θ))2fθ(x)dx dla zmiennej ciągłej;
P
x(x − µ(θ))2fθ(x) dla zmiennej dyskretnej.
W modelu bayesowskim, µ(θ) i σ2(θ) można interpretować jako warunkową wartość oczekiwaną i warunkową wariancję:
µ(θ) = E(Xi|Θ = θ), σ2(θ) = Var(Xi|Θ = θ).
Zauważmy, że µ(Θ) = E(Xi|Θ) i σ2(Θ) = Var(Xi|Θ) są zmiennymi losowymi.
Dalszą dyskusję ułatwią nam następujące oznaczenia:
m = Z
µ(θ)π(θ)dθ, s2 =
Z
σ2(θ)π(θ)dθ, a2 =
Z
(µ(θ) − m)2π(θ)dθ.
Całki w powyższych wzorach intertetujemy jako wartości oczekiwane wzglę-dem rozkładu a priori. W zasadzie trzymamy się symboli tradycyjnie uży-wanych w teorii wiarogodności. Może tylko oznaczenie a2 nie jest całkiem ortodoksyjne. Reasumując,
µ(Θ) = E(Xi|Θ), σ2(Θ) = Var(Xi|Θ), m = Eµ(Θ),
s2 = Eσ2(Θ), a2 = Varµ(Θ).
Na mocy dobrze znanej własności warunkowej wartości oczekiwanej,
EXi = EE(Xi|Θ) = Eµ(Θ) = m.
Równie znany „wzór na dekompozycję wariancji” daje
VarXi = VarE(Xi|Θ) + EVar(Xi|Θ)
= Varµ(Θ) + Eσ2(Θ)
= a2+ s2.
Dla i 6= k mamy Covθ(Xi, Xk) = 0. Istotnie, zmienne Xi i Xk są warunkowo niezależne, a więc są warunkowo nieskorelowane, dla danego Θ = θ. Stąd wynika, że
Cov(Xi, Xk) = Cov (E(Xi|Θ), E(Xi|Θ)) + ECov(Xi, Xk|Θ)
= Cov(µ(Θ), µ(Θ)) + 0
= a2.
Zastosowaliśmy „dekompozycję kowariancji”, która wygląda zupełnie podob-nie jak dla wariancji. Rówpodob-nie łatwo sprawdzić, że Cov (Xi, µ(Θ)) = a2. Otrzymane powyżej wyniki można prosto zapisać używając symbolu Iik zde-finiowanego w taki sposób:
Iik =
(1 jeśli i = k;
0 jeśli i 6= k.
4.4.3 Stwierdzenie. W modelu określonym przez 4.3.3, zmienne X1, . . . , Xn mają tę samą wartość oczekiwaną m i następującą strukturę kowariancji:
Cov(Xi, Xk) = a2+ Iiks2,
Ponadto, Eµ(Θ) = m, Varµ(Θ) = a2 i Cov (Xi, µ(Θ)) = a2.
Stwierdzenie 4.4.3 opisuje strukturę wartości oczekiwanych i kowariancji, a więc zawiera wszystko co potrzeba do rozwiązania zadania predykcji liniowej.
Przypuśćmy, że interesuje nas predykcja zmiennej losowej Xn+1, która repre-zentuje przyszłe szkody. Jeśli włączymy tę nową zmienną do naszego mo-delu, to łatwo widać, że BLP(Xn+1) = BLP(µ(Θ)). Zmienna losowa µ(Θ), powtórzmy, reprezentuje średnią wysokość szkód dla rozpatrywanego klienta.
Żeby wyznaczć BLP(µ(Θ)), musimy rozwiązać układ równań
n
X
i=1
ck(Iiks2+ a2) = a2, (k = 1, . . . , n),
c0 = m −
n
X
i=1
cim.
Przepiszmy k-te równanie w postaci
s2ck+ a2c• = a2, gdzie c•=Pn
i=1ci. Jeśli zsumujemy te równania względem k to otrzymamy s2c•+ a2nc• = a2n,
skąd
c• = a2n a2n + s2.
Niech z = c•. Teraz już łatwo wyliczyć, że ck = z/n dla k = 1, . . . , n.
Wreszcie, c0 = m − c•m = (1 − z)m.
W ten sposób wyprowadziliśmy następujący wynik:
4.4.4 TWIERDZENIE (Najlepszy liniowy predyktor). W modelu 4.3.3 mamy
BLP(µ(Θ)) = z ¯X + (1 − z)m, gdzie
X =¯ 1 n
n
X
i=1
Xi, z = a2n a2n + s2.
Liczba z nazywa się współczynnikiem zaufania (lub wiarogodności).
W Przykładzie 4.3.10 spotkaliśmy już ten sam predyktor. Ogólniej, jeśli najlepszy predyktor jest liniową funkcją obserwacji, to siłą rzeczy pokrywa się z najlepszym liniowym predyktorem, a więc musi mieć postać podaną powyżej. Modele bayesowskie, w których tak właśnie jest, to znaczy BP = BLP, nazywają się „dokładnymi modelami zaufania” (exact credibility). We wszystkich przykładach rozpatrzonych poprzednio mieliśmy do czynienia z takimi właśnie modelami.
Model Bühlmanna - Strauba.
Rozpatrzymy p jednorodnych grup danych. Dla ustalenia uwagi powiedzmy, że Xji jest średnią szkód i-tego klienta z j-tej grupy. Nasze dane mają więc postać tablicy zmiennych losowych. Zakładamy, że z każdą grupą klientów, czyli z każdym wierszem tablicy związana jest inna zmienna strukturalna:
Θ1; X11, . . . X1i, . . . X1n1, ... ... . .. ... . ..
Θj; Xj1, . . . Xji, . . . Xjnj, ... ... . .. ... . ..
Θp; Xp1, . . . Xpi, . . . Xpnp.
4.4.5 Założenie. Rozważany układ zmiennych losowych Θj; Xji (j = 1, . . . , p; i = 1, . . . , nj) ma łączną gęstość prawdopodobieństwa
f (θj); (xji) =Y
j
π(θj)Y
i
fθj(xji).
Zmienne Xji mają skończoną wartość oczekiwaną i wariancję. Innymi słowy,
• Θ1, . . . , Θp są niezależnymi zmiennymi losowymi.
• π jest gęstością prawdopodobieństwa każdej ze zmiennych Θj;
• Dla ustalonych wartości Θj = θj,
Xji są warunkowo niezależnymi zmiennymi losowymi;
fθj jest warunkową gęstością każdej ze zmiennych Xji.
Zauważmy, że warunkowy rozkład zmiennych Xji z j-tej grupy zależy tylko od Θj. Układ zmiennych losowych spełniających Założenie 4.4.5 będziemy nazywać modelem Bühlmanna - Strauba.
Funkcje µ, σ2 i liczbę m definiujemy tak jak poprzednio w terminach gęstości fθ i rozkładu a priori π. Mamy teraz
µ(Θj) = E(Xji|Θj), σ2(Θj) = Var(Xji|Θj).
Oczywiście,
m = EXji = Eµ(Θj), s2 = Eσ2(Θj), a2 = Varµ(Θj).
Opiszemy strukturę kowariancji zmiennych występujących w modelu Bühl-manna - Strauba. Wobec tego, że wiersze są niezależne, z łatwością otrzy-mujemy następujący wzór:
Cov(Xji, Xj0i0) = Ijj0(Iii0s2 + a2).
Podobnie, Cov(Xji, µ(Θj0)) = Ijj0a2.
Najlepsza predykcja liniowa. Interesuje nas głównie predykcja zmiennej losowej µ(Θj) dla pewnego ustalonego j = j0. Jest to średnia wysokość szkód w j0-tej grupie klientów. Dla ubezpieczyciela predykcja tej zmiennej ma zasadnicze znaczenie. Jest równoważna predykcji przyszłych szkód każdego z klientów należących do tej grupy (lub nawet nowego klienta, jeśli mamy podstawy zakwalifikować go właśnie do j0-tej grupy).
Zaczniemy od spostrzeżenia tyleż oczywistego, co zaskakującego. Najlepszy liniowy predyktor BLP(µ(Θj0)) obliczony na podstawie całej tablicy danych (Xji) zależy tylko od zmiennych Xj0i z j0-tej grupy. To wynika z faktu, że wiersze: (Θj; Xji, i = 1, . . . , nj) są niezależne dla j = 1, . . . , p. Stosując Twierdzenie 4.4.4 do j-tego wiersza dostajemy
BLP(µ(Θj)) = zjX¯j + (1 − zj)m, gdzie
zj = a2nj a2nj + s2
jest współczynnikiem zaufania (zależnym od grupy j), X¯j = 1
nj X
i
Xji
jest średnią danych w j-tej grupie. Oczywiście, można również odwołać się explicite do Twierdzenia 4.4.2 i wykorzystać fakt, że zmienne z różnych wier-szy są nieskorelowane, aby otrzymać w drodze czysto rachunkowej ten sam wynik.
Nasuwa się wobec tego pytanie: po co nam model uwzględniający dane dla wielu klientów, skoro najlepszy liniowy predyktor korzysta tylko z danych dotyczących jednego klienta? Odpowiedź jest bardzo prosta. Nasz predyktor zakłada znajomość trzech parametrów: m, s2 i a2. W praktyce te parametry są nieznane i musimy je estymować na podstawie danych. Dwa z nich, m i a2 opisują własności populacji klientów i mogą być estymowane tylko jeśli mamy próbkę z tej populacji! Rolę tej próbki odgrywają wiersze tablicy (Xji) – i dlatego są potrzebne. Model Bühlmanna - Strauba reprezentuje w tym sensie empiryczne podejście bayesowskie.
Istnieje pewna „hierarchia predyktorów”, której warto się przyjrzeć na przy-kładzie modelu Bühlmanna - Strauba. Załóżmy, że celem jest predykcja µ(Θj) dla ustalonego j = j0.
• BP (Best Predictor ), najlepszy predyktor: obliczenie BP wymaga zna-jomości wiarogodności fθ i rozkładu a priori π;
• BLP (Best Linear Predictor ), najlepszy liniowy predyktor: obliczenie BLP wymaga znajomości parametrów m, s2, i a2;
• BLUP (Best Linear Unbiased Predictor ), najlepszy nieobciążony li-niowy predyktor: obliczenie BLUP wymaga znajomości parametrów s2, i a2, czyli komponentów wariancyjnych; średnią globalną m esty-mujemy z danych;
• EBLUP (Empirical Best Linear Unbiased Predictor ), empiryczna wer-sja najlepszego nieobciążonego liniowego predyktora: EBLUP można obliczyć na podstawie danych; z danych estymujemy parametry m, s2, i a2
Poprzestaniemy na podaniu podstawowych wzorów na predyktory i estyma-tory używane w modelu Bühlmanna - Strauba. Wyjaśnimy sens tych wzo-rów na poziomie intuicyjnym. Pominiemy formalne definicje takich pojęć jak BLUP i BLUE (Best Linear Unbiased Estimator ), najlepszy liniowy nie-obciążony estymator oraz wyjaśnienie związku między nimi. To są tematy należące do teorii mieszanych modeli liniowych i tylko w tym kontekście można je właściwie zrozumieć.
Oznaczenia. Niech µj = µ(Θj). Jak zwykle, pomijamy wskaźniki sumowa-nia, dostatecznie jasno wynikacjące z kontekstu. Stosujemy konwencję wpro-wadzoną już poprzednio: sumowanie względem pewnego indeksu oznaczmy
„wykropkowując” ten indeks. Na przykład:
n•=X
Będziemy potrzebowali kilku skrótów używanych w analizie wariancji (ANOVA).
Przypomnijmy oznaczenia wprowadzone w Podrozdziale 4.2:
X =¯ X
Predyktory i estymatory Najlepszy liniowy predyktor BLP(µj), został obliczony przy założeniu, że znana jest globalna średnia m i zależy tylko od zmiennych Xji z j-tej grupy. Jeśli m jest nieznane, to możemy po prostu zastępujemy ten parametr odpowiednim estymatorem. Do estymacji m wy-korzystujemy dane ze wszystkich grup. Można pokazać, że najlepszy liniowy nieobciążony estymator parametru m jest dany następującym wzorem:
BLUE(m) = ˆm =X
j
zj z•
X¯j.
Interesujące jest to, że BLUE(m) jest średnią ważoną w której wagi są zwią-zane ze współczynnikami wiarogodności zj. Na pierwszy rzut oka wydawać by się mogło, że bardziej naturalne jest użycie „zwykłej” średniej ¯X. Na ogół
jednak ¯X 6= ˆm. Obie średnie są estymatorami nieobciążonymi, ale ˆm ma mniejszą wariancję (w istocie, najmniejszą spośród wszystkich estymatorów liniowych nieobciążonych). To wynika z faktu, że średnie grupowe ¯Xj mają wariancje odwrotnie proporcjonalne do zj, mianowicie Var ¯Xj = a2 = s2/nj. Wstawiając estymator BLUE(m) = ˆm w miejsce nieznanego parametru m we wzorze na BLP, otrzymujemy BLUP (możemy to prowizorycznie uznać za definicję BLUPa):
BLUP(µj) = zjX¯j + (1 − zj) ˆm.
Przejdźmy teraz do estymacji komponentów wariancyjnych s2 i a2. To jest trudniejsze zadanie i wyniki teoretyczne są mniej zadowalające, niż w przy-padku estymacji parametru m. Z pierwszym komponentem, s2 jest jeszcze nie tak trudno. Estymator
ˆ
s2 = MSB = SSW n•− p
jest nieobciążony i intuicyjnie przekonujący. Gorzej z komponentem a2. Są używane różne estymatory i nie można definitywnie powiedzieć, które z nich są „najlepsze”. Interpretacja komponentu a2podpowiada, żeby użyć „natural-nego” estymatora MSB = SSB/(p − 1). Okazuje się jednak, że ten estymator jest obciążony. Proste ale dość żmudne rachunki prowadzą do wzoru na war-tość oczekiwaną SSB w naszym modelu i pozwalają „usunąć obciążenie”. W rezultacie otrzymujemy następujący nieobciążony estymator:
ˆ a2 =
SSB − p − 1 n•− pSSW
n•
n2•−P n2j.
Powiemy, że jest to estymator otrzymany metodą ANOVA. Jego wadą jest to, że czasami przyjmuje wartości ujemne, choć jest estymatorem nieujemnego parametru a2. Z praktycznego punktu widzenia problem jest niewielki, bo można używać estymatora max(ˆa2, 0). Jednak w ten sposób otrzymujemy oczywiście estymator obciążony, a więc rezygnujemy z ważnej teoretycznej zalety estymatora ˆa2.
Zestawienie wzorów
Współczynnik zaufania dla grupy j:
zj = a2nj a2nj + s2
Najlepszy liniowy predyktor:
BLP(µj) = zjX¯j + (1 − zj)m.
Najlepszy liniowy nieobciążony predyktor:
BLUP(µj) = zjX¯j + (1 − zj) ˆm.
Najlepszy liniowy nieobciążony estymator średniej:
BLUE(m) = ˆm =X
j
zj
z•
X¯j.
Empiryczna wersja predyktora:
EBLUP(µj) = ˆzjX¯j + (1 − ˆzj) ˆm,ˆ gdzie
ˆ
zj = ˆa2nj ˆ
a2nj+ ˆs2, m =ˆˆ X
j
ˆ zj ˆ z•
X¯j.
Model Bühlmanna-Strauba jako mieszany model liniowy
Z Założenia 4.4.5 wynika, że spełnione są założenia następującego modelu, który należy do rodziny tak zwanych mieszanych modeli liniowych (Mixed Linear Models).
4.4.6 Założenie (Model 1-kierunkowej klasyfikacji z efektami losowymi).
Zmienne losowe Xji są postaci
Xji = m + αj + εji, (j = 1, . . . , p; i = 1, . . . , nj),
gdzie m jest stałą, wszystkie zmienne losowe αj i εji są nieskorelowane, Eαj = Eεji= 0, Varαj = a2, Varεji = s2.
Aby przejść od modelu bayesowskiego do modelu liniowego, czyli pokazać, że z Założenia 4.4.5 wynika Założenie 4.4.6, wystarczy napisać
εji = Xji− µ(Θj);
αj = µ(Θj) − m.
Obliczenie wartości oczekiwanych, wariancji i kowariancji tak zdefiniowanych zmiennych nie przedstawia żadnych trudności.
Interpretacja wielkości występujących w powyższym modelu jest bardzo przej-rzysta. Każdą ze zmiennych Xji rozkładamy na sumę trzech składników (w statystycznym żargonie – efektów). Liczba m jest średnią wysokością szkód w całej populacji klientów, αj jest losowym efektem związanym z przyna-leżnością klienta do j-tej grupy, zaś εji jest „błędem losowym” zależnym od grupy i od klienta. Nasze rozważania, oczywiście, żywo przypominają to, co mówiliśmy w poprzednim rozdziale o interpretacji modelu bayesowskiego.
Nic dziwnego, mamy do czynienia po prostu z opisem tego samego zjawiska w nieco innym języku.