Kilka lat temu napisałem o tym artykuł dla moich uczniów (w języku hiszpańskim), więc mogę spróbować przepisać te wyjaśnienia tutaj. Spojrzę na IRLS (iteracyjnie przeważone najmniejsze kwadraty) przez szereg przykładów o coraz większej złożoności. W pierwszym przykładzie potrzebujemy koncepcji rodziny o skali lokalizacji. Niech będzie w pewnym sensie funkcją gęstości wyśrodkowaną na zero. Możemy zbudować rodzinę gęstości, definiując
gdzie to parametr skali, a f ( x ) = f ( x ; μ , σ ) = 1f0σ>0μf0N(μ,σ)
f(x)=f(x;μ,σ)=1σf0(x−μσ)
σ>0μjest parametrem lokalizacji. W modelu błędu pomiaru, gdzie zwykle termin błędu jest modelowany jako rozkład normalny, możemy zamiast tego rozkładu normalnego użyć rodziny o skali lokalizacji skonstruowanej powyżej. Gdy jest standardowym rozkładem normalnym, powyższa konstrukcja daje rodzinę .
f0N(μ,σ)
Teraz użyjemy IRLS na kilku prostych przykładach. Najpierw znajdziemy estymatory ML (maksymalne prawdopodobieństwo) w modelu
o gęstości
Cauchy rozkład rodziny lokalizacji (więc jest to rodzina lokalizacji). Ale najpierw jakiś zapis. Estymator najmniejszych kwadratów ważony jest wyrażony przez
gdzie to niektóre wagi. Zobaczymy, że estymator ML można wyrazić w tej samej formie, za pomocąf ( y ) = 1
Y1,Y2,…,Yni.i.d
f(y)=1π11+(y−μ)2,y∈R,
μμμ∗=∑ni=1wiyi∑ni=1wi.
wiμwijakaś funkcja residuals
Funkcja prawdopodobieństwa jest podana przez
a funkcja loglikelihood jest dana przez
Jego pochodną względem jest
gdzie . pisać
ϵi=yi−μ^.
L(y;μ)=(1π)n∏i=1n11+(yi−μ)2
l(y)=−nlog(π)−∑i=1nlog(1+(yi−μ)2).
μ∂l(y)∂μ===0−∑∂∂μlog(1+(yi−μ)2)−∑2(yi−μ)1+(yi−μ)2⋅(−1)∑2ϵi1+ϵ2i
ϵi=yi−μf0(ϵ)=1π11+ϵ2 i , otrzymujemy
Znajdujemy
którym użyliśmy definicji
f′0(ϵ)=1π−1⋅2ϵ(1+ϵ2)2f′0(ϵ)f0(ϵ)=−1⋅2ϵ(1+ϵ2)211+ϵ2=−2ϵ1+ϵ2.
∂l(y)∂μ===−∑f′0(ϵi)f0(ϵi)−∑f′0(ϵi)f0(ϵi)⋅(−1ϵi)⋅(−ϵi)∑wiϵi
wi=f′0(ϵi)f0(ϵi)⋅(−1ϵi)=−2ϵi1+ϵ2i⋅(−1ϵi)=21+ϵ2i.
Pamiętając, że
otrzymujemy równanie
które jest równaniem szacunkowym IRLS. Zauważ, że
ϵi=yi−μ∑wiyi=μ∑wi,
- Wagi są zawsze dodatnie.wi
- Jeśli reszta jest duża, przypisujemy mniejszą wagę do odpowiedniej obserwacji.
Aby obliczyć estymator ML w praktyce, potrzebujemy wartości początkowej , moglibyśmy na przykład użyć mediany. Za pomocą tej wartości obliczamy resztki
i wagi
Nowa wartość jest podana przez
Kontynuując w ten sposób, definiujemy
i
Szacowana wartość na przejściu algorytmu staje się
μ^(0)
ϵ(0)i=yi−μ^(0)
w(0)i=21+ϵ(0)i.
μ^μ^(1)=∑w(0)iyi∑w(0)i.
ϵ(j)i=yi−μ^(j)
w(j)i=21+ϵ(j)i.
j+1μ^(j+1)=∑w(j)iyi∑w(j)i.
Kontynuując, aż sekwencja
zbiegnie się.
μ^(0),μ^(1),…,μ^(j),…
Teraz badamy ten proces z bardziej ogólną rodziną lokalizacji i skali, , z mniejszą ilością szczegółów. Niech będą niezależne od powyższej gęstości. Zdefiniuj także . Funkcja loglikelihood to
Pisząc , zwróć uwagę, że
i
Obliczanie pochodnej logarytmu
f(y)=1σf0(y−μσ)Y1,Y2,…,Ynϵi=yi−μσ
l(y)=−n2log(σ2)+∑log(f0(yi−μσ)).
ν=σ2∂ϵi∂μ=−1σ
∂ϵi∂ν=(yi−μ)(1ν−−√)′=(yi−μ)⋅−12σ3.
∂l(y)∂μ=∑f′0(ϵi)f0(ϵi)⋅∂ϵi∂μ=∑f′0(ϵi)f0(ϵi)⋅(−1σ)=−1σ∑f′o(ϵi)f0(ϵi)⋅(−1ϵi)(−ϵi)=1σ∑wiϵi
i tego do zera daje to samo równanie szacunkowe jak w pierwszym przykładzie. Następnie wyszukaj estymator dla :
σ2∂l(y)∂ν=====−n21ν+∑f′0(ϵi)f0(ϵi)⋅∂ϵi∂ν−n21ν+∑f′0(ϵi)f0(ϵi)⋅(−(yi−μ)2σ3)−n21ν−121σ2∑f′0(ϵi)f0(ϵi)⋅ϵi−n21ν−121ν∑f′0(ϵi)f0(ϵi)⋅(−1ϵi)(−ϵi)⋅ϵi−n21ν+121ν∑wiϵ2i=!0.
prowadząc do estymatora
W tym przypadku można również zastosować powyższy algorytm iteracyjny.
σ2^=1n∑wi(yi−μ^)2.
Poniżej podajemy przykład numeryczny z wykorzystaniem R, dla modelu podwójnego wykładniczego (o znanej skali) i danych y <- c(-5,-1,0,1,5)
. Dla tych danych prawdziwa wartość estymatora ML wynosi 0. Wartość początkowa będzie wynosić mu <- 0.5
. Jeden przebieg algorytmu to
iterest <- function(y, mu) {
w <- 1/abs(y-mu)
weighted.mean(y,w)
}
za pomocą tej funkcji możesz eksperymentować z wykonywaniem iteracji „ręcznie”. Następnie można wykonać algorytm iteracyjny
mu_0 <- 0.5
repeat {mu <- iterest(y,mu_0)
if (abs(mu_0 - mu) < 0.000001) break
mu_0 <- mu }
Ćwiczenie: Jeśli model jest rozkładem z parametrem skali pokaż iteracje według wagi
Ćwiczenie: jeśli gęstość jest logistyczna, pokaż wagi podane przez
tkσw(ϵ)=1-eϵ
wi=k+1k+ϵ2i.
w(ϵ)=1−eϵ1+eϵ⋅−1ϵ.
Na razie zostawię to tutaj, będę kontynuować ten post.