Moments of Bayesian product

We are interested in the first two moments of the product between an exponential-family likelihood with a Normal distribution.

An univariate, shape-scale exponential-family distribution is defined here as

\begin{align} f(y~{}|~{}\theta,\phi)=\exp\{(y\theta-b(\theta))/a(\phi)+c(y,\phi)\}\notag \\ \end{align}

where \(\theta\) is the mandatory shape-parameter and \(\phi\) is the optional scale-parameter. We will refer it simply as exponential-family distribution henceforth. It shows the properties

\begin{align} \mathrm{E}[y] & =\mu=b^{\prime}(\theta)\notag \\ \mathrm{Var}[y] & =\sigma^{2}=b^{\prime\prime}(\theta)a(\phi)\notag \\ \end{align}

A Normal distribution can be written in the exponential-family way as

\begin{align} \mathcal{N}(x~{}|~{}\theta,\phi)=\exp\{(x\theta-\theta^{2}/2)/\phi-(x^{2}/\phi+\log(2\pi\phi))/2\}\notag \\ \end{align}

Given any exponential-family distribution \(p(y~{}|~{}\theta,\phi)\), we want to compute

\begin{align} \int_{-\infty}^{\infty}\theta^{m}p(y~{}|~{}\theta,\phi)\mathcal{N}(\theta~{}|~{}\mu,\sigma^{2})\mathrm{d}\theta\Big{/}Z,\notag \\ \end{align}

where \(Z=\int_{-\infty}^{\infty}p(y~{}|~{}\theta,\phi)\mathcal{N}(\theta~{}|~{}\mu,\sigma^{2})\mathrm{d}\theta\) normalizes it and \(m\in\{1,2\}\), in less than a millisecond. A method that solves this problem must also provide error guarantees.

Let us thus define and derive

\begin{align} g(\theta) & =\log\left(p(y~{}|~{}\theta,\phi)\mathcal{N}(\theta~{}|~{}\mu,\sigma^{2})\right)-c(y,\phi)+\mu^{2}/(2\sigma^{2})+\log(2\pi\sigma^{2})/2\notag \\ & =\theta(y/a(\phi)+\mu/\sigma^{2})-\theta^{2}/(2\sigma^{2})-b(\theta)/a(\phi)\notag \\ g^{\prime}(\theta) & =y/a(\phi)+\mu/\sigma^{2}-\theta/\sigma^{2}-b^{\prime}(\theta)/a(\phi)\notag \\ g^{\prime\prime}(\theta) & =-\sigma^{-2}-b^{\prime\prime}(\theta)/a(\phi)\notag \\ \end{align}

Its second-order Taylor expansion, after some simplification, is

\begin{align} \tilde{g}(\theta)= & -b(\theta_{0})/a(\phi)+(b^{\prime}(\theta_{0})/a(\phi))\theta_{0}-(b^{\prime\prime}(\theta_{0})/a(\phi))\theta_{0}^{2}/2\notag \\ & +\theta(y/a(\phi)+\eta-b^{\prime}(\theta_{0})/a(\phi))+(b^{\prime\prime}(\theta_{0})/a(\phi))\theta_{0})\notag \\ & +\theta^{2}(-\tau-b^{\prime\prime}(\theta_{0}))/2\notag \\ \end{align}

Algorithm