the parameter estimation procedure can be easily visualized. Assume that given latent variables $Z_{1},\ldots,Z_{n}\overset{iid}{\sim}N(0,1)$, the data are independently generated as $X_{i}|\{Z_{1}=z_{1},\ldots,Z_{n}=z_{n}\}\sim\mathsf{Gamma}(10\cdot\sigma(a+bz_{i}))$, where $\sigma(x)=1/(1+\exp(-x))$ is the sigmoid function, and $\mathsf{Gamma}(s)$ stands for a gamma distribution with shape parameter $s$. The target is to estimate the unknown parameters $\theta=(a,b)$ from the observed data $X_{1},\ldots,X_{n}$. In our simulation, the true parameters are set to $a=2$ and $b=0.5$, and a sample size $n=100$ is used to simulate $X_{i}$.

For a single variable pair $(x,z)$, it is easy to show that the complete log-likelihood function is $L(\theta;x,z)=-z^{2}/2+(s-1)\log(x)-\log\{\Gamma(s)\}+C$, where $s=10\cdot\sigma(a+bz)$, $\Gamma(\cdot)$ is the gamma function, and $C$ is a constant. The EM algorithm is then used to solve this problem as follows. In the $k$-th M-step, we fix parameter estimate at $\theta_{k}=(a_{k},b_{k})$, and then optimize the objective function $Q(\theta;\theta_{k})=\mathbb{E}_{Z|X=x,\theta_{k}}[L(\theta;x,Z)]$ using SAGD. The next $\theta$ value is set to the optimum of $Q(\theta;\theta_{k})$.

Since for this model we can evaluate the true derivatives of $Q(\theta;\theta_{k})$ using numerical integration, it is of interest to compare SAGD with the exact gradient descent (GD) method. We set the initial value to be $\theta_{0}=(0,1)$, and run both SAGD and exact GD for $T=100$ iterations in each M-step, with a constant step 