## 3 Results
### Riemann Approximation EM: Two Applications  
#### 3.1.1 Application to a Gaussian Model with the Beta Prior  
We demonstrate the interest of the method on a example with a continuous bounded random variable following a Beta distribution \(z{\sim}Beta(\alpha,1)\), and an observed random variable following \(x{\sim}\mathcal{N}(\lambda z,\sigma^{2})\). In other words, with \(\epsilon{\sim}\mathcal{N}(0,1)\) independent of \(z\):  
\[x=\lambda z+\sigma\epsilon\,.\]  
This results in a likelihood belonging to the exponential family:  
\[h(z;\theta)=\frac{\alpha z^{\alpha-1}}{\sqrt{2\pi\sigma^{2}}}exp\left(-\frac {(x-\lambda z)^{2}}{2\sigma^{2}}\right)\,.\]  
Since \(z\) is bounded, and everything is continuous in the parameter \((\alpha,\lambda,\sigma^{2})\), this model easily verifies each of the conditions \(M1\)-\(3\). The E step with this model involves the integral \(\int_{z}z^{\alpha}exp\left(-\frac{(x-\lambda z)^{2}}{2\sigma^{2}}\right)dz\), a fractional moment of the Gaussian distribution. Theoretical formulas exists for these moments, see [40], however they involve [PERSON]'s confluent hypergeometric functions, which are infinite series. Instead, we use the Riemann approximation to run the EM algorithm with this model: \(\tilde{h}_{n}(z;\theta):=h(\lfloor\varphi(n)z\rfloor/\varphi(n);\theta)\). As done previously, we take, without loss of generality, \(\varphi(n):=n\) for the sake of simplicity. The E step only involves the \(n\) different values taken by the step function probabilities \(h(\lfloor nz\rfloor/n;\theta)\):  
\[\tilde{p}_{\theta,n}^{(i)}\left(\frac{k}{n}\right)=\frac{h^{(i)}(\frac{k}{n}; \theta)}{\frac{1}{n}\sum_{l=0}^{n-1}h^{(i)}(\frac{l}{n};\theta)}\,.\]  
where the exponent \((i)\) indicates the index of the observation \(x^{(i)}\). To express the corresponding M step in a digest way, let us define the operator \(\Psi^{(i)}:\mathbb{R}^{[0,1]}\longrightarrow\mathbb{R}\) such that, for any \(f:[0,1]\longrightarrow\mathbb{R}\):  
\[\Psi^{(i)}\circ f=\sum_{k=0}^{n-1}\tilde{p}_{\theta,n}^{(i)}\left(\frac{k}{n} \right)\int_{z=k/n}^{(k+1)/n}f(z)dz\,.\]  
Then, the M step can be expressed as:  
\[\frac{1}{\hat{\alpha}}=-\frac{1}{N}\sum_{i=1}^{N}\Psi^{(i)}\circ \ln(z)\,,\] \[\hat{\lambda}=\frac{\sum_{i=1}^{N}\Psi^{(i)}\circ(x^{(i)}z)}{\sum _{i=1}^{N}\Psi^{(i)}\circ z^{2}}\,, \tag{13}\] \[\tilde{\sigma}^{2}=\frac{1}{N}\sum_{i=1}^{N}\Psi^{(i)}\circ\left( x^{(i)}-\hat{\lambda}z\right)^{2}.\]  
where we took the liberty of replacing \(f\) by \(f(z)\) in these equations for the sake of simplicity. Here \(N\) is the total number of observations: \(x:=(x^{(1)},...,x^{(N)})\) idl. We test this algorithm on synthetic data. With real values \(\alpha=2,\lambda=5,\sigma=1.5\), we generate a dataset with \(N=100\) observations and run the Riemann EM with random initialisation. This simulation is ran 2000 times. We observe that the Riemann EM is indeed able to increase the likelihood, despite the EM being originally intractable. On Figure 1, we display the average trajectory, with standard deviation, of the negative log-likelihood \(-\ln(g(\theta))\) during the Riemann EM procedure. The profile is indeed decreasing. The standard deviation around the average value is fairly high, since each run involves a different dataset and a different random initialisation, hence different value of the likelihood, but the decreasing trend is the same for all of the runs. We also display the average relative square errors on the parameters at the end of the algorithm. They are all small, with reasonablysmall standard deviation, which indicates that the algorithm consistently recovers correctly the parameters.