## 3 Estimation Procedure Based on StEM
### The Estimation Procedure
[label=Step 1]
2. simulate \(\mathbf{y}_{cen,i}^{*}\) from \(\text{TN}(d,\boldsymbol{\mu}_{i}^{(k-1)},\Sigma_{i}^{(k-1)})\), a multivariate truncated normal distribution with * \(d\) the lower bound of the truncation, * mean \(\boldsymbol{\mu}_{i}^{(k-1)}=(\mu_{i,1}^{(k-1)},\ldots,\mu_{i,n_{i}}^{(k-1)})\) where \(\mu_{ij}^{(k-1)}=g((t_{ij}-\tau_{i}^{(k-1)})^{-},\mathbf{a}_{i}^{(k-1)})+h((t_ {ij}-\tau_{i}^{(k-1)})^{+},\mathbf{b}_{i}^{(k-1)})\), \(j=1,\ldots,n_{i}\), * and variance \(\Sigma_{i}^{(k-1)}=\sigma_{e}^{(k-1)}I_{n_{i}\times n_{i}}\), where \(I_{n_{i}\times n_{i}}\) is a \(n_{i}\) by \(n_{i}\) identity matrix, and, independently, sample \(\eta\) from the uniform (0, 1) distribution;
3. calculate \(\rho=\frac{t_{i}^{(0)}\left(\boldsymbol{\theta}^{(k-1)}\left|\left(\mathbf{y }_{cen,i}^{(k-1)},\mathbf{y}_{cen,i}^{(k-1)},\mathbf{b}_{i}^{(k-1)},\mathbf{ y}_{cen,i}^{(k-1)}\right)\right.\right)}{t_{i}^{(0)}\left(\boldsymbol{\theta}^{(k-1)} \left|\left(\mathbf{y}_{cen,i}^{(k-1)},\mathbf{y}_{cen,i}^{(k-1)},\mathbf{a}_ {i}^{(k-1)}\right)\right.\right.\);
4. if \(\eta\leq\rho\), we update \(\mathbf{y}_{cen,i}^{(k)}\) by \(\mathbf{y}_{cen,i}^{*}\); otherwise, \(\mathbf{y}_{cen,i}^{(k)}=\mathbf{y}_{cen,i}^{(k-1)}\). The maximization step of the StEM algorithm involves maximizing the log-likelihood \(\sum_{i=1}^{n}l_{c}^{(i)}(\boldsymbol{\theta}|\mathbf{y}_{obs,\boldsymbol{y }},\boldsymbol{y}_{cen,i},\mathbf{a}_{i},\mathbf{b}_{i},\tau_{i})\) to update parameters \(\{\boldsymbol{\alpha},\boldsymbol{\beta},\mathbf{A},\mathbf{B},\sigma_{e}^{2},\pi,\lambda\}\) under the current imputed missing data. Since \((\mathbf{y}_{cen,i},\mathbf{a}_{i},\mathbf{b}_{i},\tau_{i})\) are regarded as data, the complete log-likelihood no longer involves integrals, which substantially simplifies the maximization. Solving the corresponding score functions yields the following estimations:  
\[\boldsymbol{\alpha} =\frac{1}{n}\sum_{i=1}^{n}\mathbf{a}_{i},\quad\boldsymbol{\beta }=\frac{1}{n}\sum_{i=1}^{n}\mathbf{b}_{i},\] \[\mathbf{A} =\frac{1}{n}\sum_{i=1}^{n}(\mathbf{a}_{i}-\boldsymbol{\alpha})^{ T}(\mathbf{a}_{i}-\boldsymbol{\alpha}),\quad\mathbf{B}=\frac{1}{n}\sum_{i=1}^{n}( \mathbf{b}_{i}-\boldsymbol{\beta})^{T}(\mathbf{b}_{i}-\boldsymbol{\beta}),\] \[\sigma_{e}^{2} =\frac{1}{n}\sum_{i=1}^{n}\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}\{y_{ ij}-g((t_{ij}-\tau_{i})^{-},\mathbf{a}_{i})-h((t_{ij}-\tau_{i})^{+},\mathbf{b}_{i})\}^{2}.\]  
The likelihood function represented by the joint probability density and mass function of the ZIE distribution can be written as \(L(\tau_{1},\tau_{2},\ldots,\tau_{n};\pi,\lambda)=\prod_{\tau_{i}=0}\pi\prod_{ \tau_{j}\geq 0}\)\((1-\pi)\lambda\exp(-\lambda\tau_{i})\). Denote \(I=\left\{ \begin{array}{ll}1,&\text{if }\tau_{i}=0\\ 0,&\text{other}\end{array} \right.\), we have \(\log L=\sum_{i=1}^{n}I\log(\pi)+\sum_{i=1}^{n}(1-I)\log[(1-\pi)\lambda\exp(- \lambda\tau_{i})]\). Solving the score function \(\partial\log L/\partial\pi=0\) and \(\partial\log L/\partial\lambda=0\), yields  
\[\pi=\frac{1}{n}\sum_{i=1}^{n}I,\quad\lambda=\Big{(}1-\sum_{i=1}^{n}I\Big{)} \Big{/}\Big{(}\sum_{i=1}^{n}\tau_{i}\Big{)}.\]  
In general, increasing the number of random effects does not substantially increase the complexity of the maximization step. However, the imputation step will be more complicated, as this increases the dimensions of missing data that need to be imputed. As with the likelihood defined in (3), the Fisher information matrix of the ZIE-based nonlinear random change point model has no closed-form solution. To obtain the variance-covariance matrix of the MLE \(\hat{\mathbf{\theta}}\), we consider the following approximate formula in [23]. Denote the score function of the complete-data likelihood by \(S_{c}^{(i)}=\partial l_{c}^{(i)}/\partial\mathbf{\theta}\).