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

Model regresji Poissona napompowany zerem jest definiowany dla próbki przez i zakłada ponadto, że parametry $(y_1,\ldots,y_n)$

Y_{ja} = {\begin{cases} 0 & z prawdopodobieństwem p_{ja} + (1 - p_{ja}) {mi}^{- λ_{ja}} \\ k & z prawdopodobieństwem (1 - p_{ja}) {mi}^{- λ_{ja}} λ_{ja}^{k} / k! \end{cases}

$Y_i = \begin{cases} 0 & \text{with probability} \ p_i+(1-p_i)e^{-\lambda_i}\\ k & \text{with probability} \ (1-p_i)e^{-\lambda_i} \lambda_{i}^{k}/k! \end{cases}$

zaspokoić

λ = (λ_{1}, \dots, λ_{n})

$\mathbf{\lambda} = (\lambda_1, \dots, \lambda_n)$

p = (p_{1}, \dots, p_{n})

$\textbf{p} = (p_1, \dots, p_n)$

\begin{aligned} \log (λ) & = b β \\ logit (p) & = \log (p / (1 - p)) = sol γ . \end{aligned}

$\eqalign{ \log(\mathbf{\lambda}) &= \textbf{B} \beta \\ \text{logit}(\textbf{p}) &= \log(\textbf{p}/(1-\textbf{p})) = \textbf{G} \mathbf{\gamma}. }$

\begin{aligned} L. (γ, β; y) & = \sum_{y_{ja} = 0} \log ({mi}^{{sol}_{ja} γ} + \exp (- {mi}^{b_{ja} β})) + \sum_{y_{ja} > 0} (y_{ja} b_{ja} β - {mi}^{b_{ja} β}) \\ - \sum_{ja = 1}^{n} \log (1 + {mi}^{{sol}_{ja} γ}) - \sum_{y_{ja} > 0} \log (y_{ja}!) \end{aligned}

$\eqalign{ L(\gamma,\mathbf{\beta}; \mathbf{y}) &= \sum_{y_i=0} \log(e^{G_i \gamma}+\exp(-e^{\textbf{B}_i \mathbf{\beta}})) +\sum_{y_i >0} (y_i \textbf{B}_i \mathbf{\beta}-e^{\textbf{B}_i \mathbf{\beta}})\\ &\quad -\sum_{i=1}^{n} \log(1+e^{G_{i} \gamma})-\sum_{y_i >0} \log(y_{i}!)}$

$\mathrm{B}$ $\mathrm{G}$

$Z_i = 1$ $Y_i$ $Z_i = 0$ $Y_i$

L. (γ, β; y, z) = \sum_{ja = 1}^{n} \log (fa (z_{ja} | γ)) + \sum_{ja = 1}^{n} \log (fa (y_{ja} | z_{ja}, β))

$L(\gamma,\mathbf{\beta}; \mathbf{y}, \mathbf{z}) = \sum_{i=1}^{n} \log(f(z_i|\mathbf{\gamma}))+\sum_{i=1}^{n} \log(f(y_i|z_i, \mathbf{\beta}))$

= \sum_{ja = 1}^{n} z_{ja} ({sol}_{ja} γ - \log (1 + {mi}^{{sol}_{ja} γ})) + - \sum_{ja = 1}^{n} (1 - z_{ja}) \log (1 + {mi}^{{sol}_{ja} γ}) + \sum_{ja = 1}^{n} (1 - z_{ja}) [y_{ja} b_{ja} β - {mi}^{b_{ja} β} - \log (y_{ja}!)]

$= \sum_{i=1}^{n} z_{i} (\textbf{G}_i \gamma-\log(1+e^{G_{i} \gamma}))+ -\sum_{i=1}^{n} (1-z_{i})\log(1+e^{G_{i} \gamma})+ \sum_{i=1}^{n} (1-z_i)[y_{i} \textbf{B}_i \beta-e^{\textbf{B}_i \beta} - \log(y_{i}!)]$

z_{i} = 0

$z_i=0$

z_{i} = 1

$z_i=1$

$Z_i = 0$ $Z_i = 1$

— Damien
źródło

f

$f$

f

$f$

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

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

$z_i$

$k^{th}$ $z_i$ $(k-1)^{th}$

$\lambda$ $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 $\beta$ i dlatego $\lambda_i$ , a następnie zmaksymalizować:

$\sum (\mathbb{E}z_i\log{p_i} + (1-\mathbb{E}z_i)\log{(1-p_i)})$

w odniesieniu do wektora współczynnika macierzy $\mathbf{G}$ uzyskać szacunki $p_i$ . Oczekiwane wartości $\mathbb{E}z_i = p_i/(p_i+(1-p_i)\exp{(-\lambda_i)})$ , 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ę