## 3 The Epsom-Hyb Estimator
### Estimating Truncated Gaussian Probabilities  
Despite the prevalence of Gaussian densities in statistical modeling, multivariate Gaussian probabilities remain difficult to compute, as they typically require high-dimensional integration. Numerical integration [29] is a common solution, but it is prohibitively expensive in the number of points required, thus preventing it from being a scalable solution in the high-dimensional space in which we wish to operate. A more recent alternative that specifically targets truncated normal distributions is the minimax tilting method [30], and while we used this method in the initial development of our algorithm, it failed to produce accurate results beyond trivial settings. This ultimately led us to consider more computationally efficient integration techniques, such as expectation propagation [31], which is widely used to compute approximate integrals. To better understand how the expectation propagation (EP) algorithm can be applied to the intractable Gaussian integral in the previous section, we lay out some preliminary groundwork for the EP algorithm. Starting with the Gaussian distribution \(p_{0}(x)=\mathcal{N}(x\mid m,K)\), we define the un-normalized truncated distribution \(p(x)=p_{0}(x)\mathbb{1}_{\mathcal{A}}(x)\). The Gaussian probability of interest can be written as:  
\[F(\mathcal{A})=\int_{\mathcal{A}}p_{0}(x)dx=\int p(x)dx. \tag{7}\]  
[PERSON] et al. [32] propose a framework for approximating \(F(\mathcal{A})\) that consists of two main approximations. The first approximation is a tractable distribution \(q\) that minimizes the Kullback-Leibler (KL) divergence between \(p\) and \(q\), denoted \(D(p\,||\,q)\). However, the intractability of \(p\) complicates the minimization and subsequently requires another approximation, for which we instead work with a simplified representation of \(p\). Namely, we replace \(p(x)\) with the product of a prior distribution \(p_{0}(x)\) and factors \(t_{i}(x)\), such that  
\[p(x)=p_{0}(x)\prod_{i}t_{i}(x).\]In order to facilitate tractability, we assume \(t_{i}(x)=\mathbb{1}\left\{a_{i}<x_{i}<b_{i}\right\}\), where \(x_{i}\) is the \(i\)-th coordinate of \(x\), and \(a_{i},b_{i}\) are simply the lower and upper bounds of integration, respectively. Then, the target integral in Equation (7) can be written as  
\[F(\mathcal{A})=\int p(x)dx=\int p_{0}(x)\prod_{i=1}^{d}t_{i}(x_{i})dx_{i}.\]  
We proceed to approximate each of the intractable factors \(t_{i}\) with a tractable, un-normalized Gaussian \(\tilde{I}_{i}(x)\), which produces the final approximation \(q\) of \(p\). More specifically, we take \(q\) to mirror the product form of \(p\),  
\[q(x)=p_{0}(x)\prod_{i}\tilde{t}_{i}(x)=p_{0}(x)\prod_{i}\tilde{Z}_{i}\, \mathcal{N}\Big{(}x\mid\tilde{\mu}_{i},\tilde{\sigma}_{i}^{2}\Big{)}=Z\mathcal{ N}(x\mid\mu,\Sigma),\]  
where the parameters of these unnormalized Gaussian distributions, \(\{\tilde{\mu}_{i},\tilde{\sigma}_{i}^{2},\tilde{Z}_{i}\}\), admit closed form updates that are the result of an iterative moment matching scheme. From this, we observe that by estimating the normalizing constant of \(q\), we also obtain an approximation for the normalizing constant of the target distribution \(p\). See Equations (21)-(23) from [PERSON] et al. (2018) for the closed-form updates for each of these parameters and more details regarding the relevant notation. Essentially, the EP algorithm iteratively constructs the approximating distribution \(q(x)\) to minimize \(D\Big{(}t_{i}q^{\setminus i}\,||\,\tilde{I}_{i}q^{\setminus i}\Big{)}\), which in turn approximately minimizes \(D(p\,||\,q)\). Here, \(q^{\setminus i}(x)=q(x)/\tilde{t}_{i}(x)\) is defined as the cavity distribution.