W jaki sposób wykorzystujesz algorytm EM do obliczania MLE dla sformułowania zmiennej utajonej modelu Poissona z napompowaniem zerowym?

10

Model regresji Poissona napompowany zerem jest definiowany dla próbki przez Y i = { 0 z prawdopodobieństwem p i + ( 1 - p i ) e - λ i k z prawdopodobieństwem ( 1 - p i ) e - λ i λ k i / k ! i zakłada ponadto, że parametry λ =(y1,,yn)

Yja={0z prawdopodobieństwem pja+(1-pja)mi-λjakz prawdopodobieństwem (1-pja)mi-λjaλjak/k!
i P = ( P 1 , ... , p n ) zaspokoićλ=(λ1,,λn)p=(p1,,pn)

log(λ)=bβlogit(p)=log(p/(1-p))=solγ.

L.(γ,β;y)=yja=0log(misoljaγ+exp(-mibjaβ))+yja>0(yjabjaβ-mibjaβ)-ja=1nlog(1+misoljaγ)-yja>0log(yja!)

bsol

Zja=1YjaZja=0Yja

L.(γ,β;y,z)=ja=1nlog(fa(zja|γ))+ja=1nlog(fa(yja|zja,β))

=ja=1nzja(soljaγ-log(1+misoljaγ))+-ja=1n(1-zja)log(1+misoljaγ)+ja=1n(1-zja)[yjabjaβ-mibjaβ-log(yja!)]
zja=0zja=1

Zja=0Zja=1

Damien
źródło
fa
fa

Odpowiedzi:

11

Źródłem trudności, na które napotykasz, jest zdanie:

Następnie za pomocą algorytmu EM możemy zmaksymalizować drugie prawdopodobieństwo dziennika.

zja

kthzja(k-1)th

λp

# Generate data
# Lambda = 1,  p(zero) = 0.1
x <- rpois(10000,1)
x[1:1000] <- 0

# Sufficient statistic for the ZIP
sum.x <- sum(x)

# (Poor) starting values for parameter estimates
phat <- 0.5
lhat <- 2.0

zhat <- rep(0,length(x))
for (i in 1:100) {
  # zhat[x>0] <- 0 always, so no need to make the assignment at every iteration
  zhat[x==0] <- phat/(phat +  (1-phat)*exp(-lhat))

  lhat <- sum.x/sum(1-zhat) # in effect, removing E(# zeroes due to z=1)
  phat <- mean(zhat)   

  cat("Iteration: ",i, "  lhat: ",lhat, "  phat: ", phat,"\n")
}

Iteration:  1   lhat:  1.443948   phat:  0.3792712 
Iteration:  2   lhat:  1.300164   phat:  0.3106252 
Iteration:  3   lhat:  1.225007   phat:  0.268331 
...
Iteration:  99   lhat:  0.9883329   phat:  0.09311933 
Iteration:  100   lhat:  0.9883194   phat:  0.09310694 

W twoim przypadku na każdym kroku wykonasz ważoną regresję Poissona, w której wagi mają 1-zhatuzyskać oszacowaniaβ i dlatego λja, a następnie zmaksymalizować:

(mizjalogpja+(1-mizja)log(1-pja))

w odniesieniu do wektora współczynnika macierzy sol uzyskać szacunki pja. Oczekiwane wartościmizja=pja/(pja+(1-pja)exp(-λja)), ponownie obliczany przy każdej iteracji.

Jeśli chcesz to zrobić dla rzeczywistych danych, w przeciwieństwie do zwykłego zrozumienia algorytmu, pakiety R już istnieją; Oto przykład http://www.ats.ucla.edu/stat/r/dae/zipoisson.htm z wykorzystaniem psclbiblioteki.

EDYCJA: Powinienem podkreślić, że to, co robimy, to maksymalizacja oczekiwanej wartości prawdopodobieństwa dziennika kompletnych danych, NIE maksymalizowanie prawdopodobieństwa dziennika kompletnych danych z podłączonymi oczekiwanymi wartościami brakujących danych / zmiennych ukrytych. Tak się składa, jeśli prawdopodobieństwo dziennika pełnych danych jest liniowe w brakujących danych, ponieważ tutaj są dwa podejścia są takie same, ale w przeciwnym razie nie są.

łucznik
źródło
@Coke, powinieneś dodać te informacje jako własną odpowiedź uzupełniającą, a nie zmieniać istniejącą odpowiedź. Ta edycja nie powinna była zostać zatwierdzona.
gung - Przywróć Monikę