Model-based Clustering of Count Processes

Abstract

A model-based clustering method based on Gaussian Cox process is proposed to address the problem of clustering of count process data. The model allows for nonparametric estimation of intensity functions of Poisson processes, while simultaneous clustering count process observations. A logistic Gaussian process transformation is imposed on the intensity functions to enforce smoothness. Maximum likelihood parameter estimation is carried out via the EM algorithm, while model selection is addressed using a cross-validated likelihood approach. The proposed model and methodology are applied to two datasets.

Appendices

Appendix A: Consistency of Penalized MLE

Sufficient conditions for $\hat {P}$ to be strongly consistent are given in Kiefer and Wolfowitz (1956) and are stated below.

• (C1) Identifiability: Let H(y_i; α_i, P) be the cumulative distribution function of h(y_i; α_i, P). If for any α_i, H(y_i; α_i, P) = H(y_i; α_i, P^∗) for all y_i, then D(P, P^∗) = 0.

• (C2) Continuity: The component parameter space $\mathcal {F}$ is a closed set. For all y_i, α_i and any given P₀,

  $$ \begin{array}{@{}rcl@{}} \lim_{P \rightarrow P_{0}} h(y_{i};{\alpha}_{i},P) = h(y_{i};{\alpha}_{i},P_{0}) {.} \end{array} $$

• (C3) Finite Kullback–Leibler Information: For any P≠P^∗, there exists an 𝜖 > 0 such that

  $$ \begin{array}{@{}rcl@{}} E^{*}[ \log(h(Y_{i}; {\alpha}_{i}, B_{\epsilon}(P) ) / h(Y_{i}; {\alpha}_{i}, P^{*}) ) ]^{+} < \infty { ,} \end{array} $$

  where B_𝜖(P) is the open ball of radius 𝜖 > 0 centered at P with respect to the metric D, and

  $$ \begin{array}{@{}rcl@{}} h(Y_{i}; {\alpha}_{i}, B_{\epsilon}(P)) = \sup_{\tilde{P} \in B_{\epsilon}(P) } h(Y_{i}; {\alpha}_{i}, \tilde{P}) { .} \end{array} $$

• (C4) Compactness: The definition of the mixture density h(y_i; α_i, P) in P can be continuously extended to a compact space $\overline { \mathbb {P} }$.

Identifiability of the mixture model is established in Proposition 3.1. In the case of mixture of Poisson processes, the component parameter space $\mathcal {F}$ is clearly closed. To compactify the space of mixing distributions $\mathbb { P}$, we notice that the distance between any two mixing distributions is bounded above by:

$$ \begin{array}{@{}rcl@{}} {\int}_{\mathcal{F} } \exp(-|\mathbf{f}|) d \mathbf{f} < \infty {.} \end{array} $$

We can extend $\mathbb {P}$ to $\overline {\mathbb { P}}$ by including all sub-distributions ρP for any ρ ∈ [0, 1) as in Chen (2017).

We now show the continuity of h(y_i; α_i, P) on $\overline { \mathbb {P} }$. Recall that $D(P_{m},P) \rightarrow 0$ if and only if $P_{m} \rightarrow P$ in distribution if only if ${\int \limits } u(\mathbf {f}) d P_{m}(\mathbf {f}) \rightarrow {\int \limits } u(\mathbf {f}) dP_{0}(\mathbf {f})$ for all bounded and continuous function u(⋅).

For any given $\mathbf {f}_{0} \in \mathcal {F}$ and α_i > 0, we clearly have $\lim _{\mathbf {f} \rightarrow \mathbf {f}_{0} } h(y_{i};{\alpha }_{i},\mathbf {f}) = h(y_{i};{\alpha }_{i},\mathbf {f}_{0})$, and $\lim _{| \mathbf {f} | \rightarrow \infty } h(y_{i}; {\alpha }_{i}, \mathbf {f} ) = 0$. Hence, h(y_i; α_i, f) is continuous and bounded on $\mathcal {F}$. Suppose that $P_{m} \rightarrow P_{0} \in \overline {\mathbb {P}}$ in distribution, we have that:

$$ \begin{array}{@{}rcl@{}} h(y_{i};{\alpha}_{i},P_{m}) = \int h(y_{i};{\alpha}_{i},\mathbf{f}) dP_{m}(\mathbf{f}) \rightarrow \int h(y_{i};{\alpha}_{i},\mathbf{f}) dP_{0}(\mathbf{f}) = h(y_{i};{\alpha}_{i},P_{0}) \end{array} $$

which shows that h(y_i; α_i, P) is continuous in $\mathbb {P}$.

To prove finite Kullback–Leibler information, we need the extra sufficient condition that all scaling factors {α_i}_i are bounded above, that is, ${\alpha }_{i} < \alpha ^{(M)} < \infty$ almost surely ∀i. Let f₀ be a support point of P^∗. There must be a positive constant δ such that:

$$ \begin{array}{@{}rcl@{}} h(y_{i};{\alpha}_{i},P^{*}) \ge \delta \left( \exp(-{\alpha}_{i} ) \prod\limits_{k=1}^{m} \frac{ ({\alpha}_{i} {\lambda}_{k})^{y_{k}} }{y_{k}!} \right) p(\mathbf{f}_{0}) { .} \end{array} $$

uniformly for all α_i < α^(M). Therefore, we have that:

$$ \begin{array}{@{}rcl@{}} E^{*}(\log h(y_{i};{\alpha}_{i},P^{*} ) )\! \ge \!\sum\limits_{k=1}^{m} \log({\alpha}_{i} {\lambda}_{k}) E^{*}(Y_{i,k}) - \! \sum\limits_{k=1}^{m} E^{*}(\log Y_{i,k}!) - {\alpha}_{i} + \log(\delta) + \log p(\mathbf{f}_{0}) { .} \end{array} $$

Since α_i < α^(M) almost surely, both E^∗(Y_{i, k}) and $E^{*}(\log (Y_{i,k}!))$ are finite for all k. Hence, $E^{*}(\log h(y_{i};{\alpha }_{i},P^{*}) ) > - \infty$. Since h(y_i; α_i, P^∗) is bounded from above, $E^{*}(\log h(y_{i};{\alpha }_{i},P^{*}) ) < \infty$. Therefore, $E^{*}|\log h(y_{i};{\alpha }_{i},P^{*}) | < \infty$.

Since $h(Y_{i};{\alpha }_{i},B_{\epsilon (P)}) < c < \infty$ for any P, 𝜖 and α_i,

$$ \begin{array}{@{}rcl@{}} E^{*}[ \log(h(Y_{i}; {\alpha}_{i}, B_{\epsilon}(P) ) / h(Y_{i};{\alpha}_{i}, P^{*}) ) ]^{+} < \log(c) - E^{*}(\log h(Y_{i};{\alpha}_{i}, P^{*})) < \infty { .} \end{array} $$

Therefore, we have shown that the four conditions above are satisfied and the penalized MLE under mixture of Poisson process is strongly consistent.

Appendix B: Derivation of EM Algorithm

Recall the complete data log-likelihood defined in Eq. 10. The E-step requires the computation of the expected value of the complete data log-likelihood function with respect to the conditional distribution of Z given x and current estimates of parameters $\hat {\theta }^{(t)}$. This conditional expression can be expressed as:

$$ \begin{array}{@{}rcl@{}} E_{Z|x,\hat{\theta}^{(t)}}(\log L(\theta;x,Z)) &=& \sum\limits_{i=1}^{n} \sum\limits_{g=1}^{G} \hat{\pi}_{i}^{g} \left[ \log(\tau_{g}) - {\alpha}_{i} + \sum\limits_{j=1}^{n_{i}} \log({\alpha}_{i}) + \log({\lambda}_{g}(x_{ij})) \right] \\ &&+ {\sum}_{g=1}^{G} \log(p(\mathbf{f}_{g})) \end{array} $$

(16)

The M-step finds the parameters θ that maximize the expression above. It is straightforward to show (by differentiation) that the values of α and τ that maximize the expression above are given by Eqs. 12 and 13 respectively.

To optimize Eq. 16 with respect to f_g, we write f_g = (f_g,1,⋯ , f_{g, m})^T and retain terms that involve f_g for g = 1,⋯ , G to obtain:

$$ \begin{array}{@{}rcl@{}} Q & \equiv & \sum\limits_{i=1}^{n} \hat{\pi}_{i}^{g} \sum\limits_{j=1}^{n_{i}} \log({\lambda}_{g}(x_{ij})) + \log(p(\mathbf{f}_{g}) ) \\ & = & \sum\limits_{i=1}^{n} {\pi_{i}^{g}} \sum\limits_{j=1}^{n_{i}} \log \left[ \frac{ \exp(f_{g}(x_{i,j}) ) }{ {\sum}_{k=1}^{m} \exp(f_{g,k}) } \right] - \frac{1}{2} \mathbf{f}_{g}^{T} K^{-1} \mathbf{f}_{g} + const \\ & = & \sum\limits_{i=1}^{n} \hat{\pi}_{i}^{g} \sum\limits_{j=1}^{n_{i}} f_{g}(x_{i,j}) - \sum\limits_{i=1}^{n} \hat{\pi}_{i}^{g} n_{i} \log\left[ \sum\limits_{k=1}^{m} \exp(f_{g,k}) \right]- \frac{1}{2} \mathbf{f}_{g}^{T} K^{-1} \mathbf{f}_{g} + const \\ & = & \mathbf{b}^{T} \mathbf{f}_{g} - \sum\limits_{i=1}^{n} \hat{\pi}_{i}^{g} n_{i} \log\left[ \sum\limits_{k=1}^{m} \exp(f_{g,k}) \right]- \frac{1}{2} \mathbf{f}_{g}^{T} K^{-1} \mathbf{f}_{g} + const \end{array} $$

where b = (b₁,⋯ , b_m)^T with y_k defined below.

$$ \begin{array}{@{}rcl@{}} b_{k} = \sum\limits_{i=1}^{n} \hat{\pi}_{i}^{g} \sum\limits_{j=1}^{n_{i}} I(x_{ij} \in [ s_{k-1}, s_{k} )) \end{array} $$

where [s_k− 1, s_k) is the k th interval of the m equally spaced intervals of [0, T).

For g = 1,⋯ , G, we now optimize Q with respect to f_g, and as analytical expression does not exist, we use Newton’s method. For each g, the Jacobian J_g(Q) of Q can be written as:

$$ \begin{array}{@{}rcl@{}} J_{g}(Q) = \mathbf{y} - \mathbf{u}_{g}\sum\limits_{i=1}^{n} \hat{\pi}_{i}^{g} n_{i} - \frac{1}{2}\mathbf{f}_{g}^{T}(K^{-1}+K^{-T}) \end{array} $$

where u_g = (u_g,1,⋯ , u_{g, m})^T with

$$ \begin{array}{@{}rcl@{}} u_{g,l} = \frac{\exp(f_{g,l})}{{\sum}_{k=1}^{m} \exp(f_{g,k})} \end{array} $$

for l = 1,⋯ , m. The Hessian H_g(Q) can be written as:

$$ \begin{array}{@{}rcl@{}} H_{g}(Q) = - (\text{diag}(\mathbf{u}_{g}) - \mathbf{u}_{g} \mathbf{u}_{g}^{T}) \sum\limits_{i=1}^{n} {\pi_{i}^{g}} n_{i} - \frac{1}{2}(K^{-1}+K^{-T}) \end{array} $$