# Recurrence formula for the moments of a half-gaussian distribution (on R+)

I am trying to compute an integral that looks like the moments of a Gaussian $$\mathcal{N}(\mu, \sigma^2)$$, but the main difference is that we only integrate over R+ and not R. I believe we could call this the moments of a Half-Gaussian distribution :

$$I_m = \int_0^{\infty}x^m \exp\left(-\frac{1}{2\sigma^2}(x-\mu)^2\right)\mathrm{d}x, \quad m\in\mathbb{N}~.$$

I was wondering if there is a (recurrence) formula of $$I_m$$ for any $$m\in \mathbb{N}$$.

So far I found that:

• $$m=0$$: $$I_0 = \dfrac{\operatorname{erf}\left(\frac{{\mu}}{\sqrt{2}\,{\sigma}}\right)+1}{2}$$
• $$m=1$$: $$I_1 = \dfrac{\sqrt{2}\operatorname{\Gamma}\left(1,\frac{{\mu}^2}{2{\sigma}^2}\right)\,{\sigma}+\left(2\sqrt{{\pi}}-\operatorname{\Gamma}\left(\frac{1}{2},\frac{{\mu}^2}{2{\sigma}^2}\right)\right){\mu}}{2\sqrt{{\pi}}}$$
• $$m=2$$: $$I_2 = \dfrac{2^\frac{3}{2}\operatorname{\Gamma}\left(1,\frac{{\mu}^2}{2{\sigma}^2}\right)\,{\mu}{\sigma}+\left(2\sqrt{{\pi}}- \operatorname{\Gamma}\left(\frac{1}{2},\frac{{\mu}^2}{2{\sigma}^2}\right) + \left(2\sqrt{{\pi}}-2\operatorname{\Gamma}\left(\frac{3}{2},\frac{{\mu}^2}{2{\sigma}^2}\right)\right){\sigma}^2\right){\mu}^2}{2\sqrt{{\pi}}}$$
• ...

with $$\Gamma$$ the incomplete Gamma function. This led me to the formula for $$m>1$$: $$I_m = T_{m,\sigma} \mu^{m-1}\sigma + T_{m,\sigma^2} \mu^{m-2}\sigma^2 + \dots + T_{m,\sigma^{m-1}} \mu\sigma^{m-1} + T_{m,\sigma^m}\sigma^m = \sum_{n=1}^m T_{m,\sigma^n} \mu^{m-n}\sigma^n$$

I managed to find a recurrence formula for $$T_{m,\sigma^n}$$ as a function of $$m$$ for $$n=1,2,3,4$$. However, I stopped there because it didn't seem to lead me to a general formula for any $$n$$. Also, the results I found are assuming that $$\mu>0$$, which is quite restrictive.

Given the nice form of the integral $$I_m$$ and its close relationship to the moments of a Gaussian, I was wondering if someone might know how to compute it and help me out...

• Your distribution is not really a Half-Gaussian if $\mu \not=0$ but a Truncated Gaussian at least if you rescale the density to integrate to $1$ Commented Mar 17 at 12:32
• The key is to use your Calculus skills to help the computer. You can relate the formula for negative $\mu$ to that for $-\mu,$ so the sign of $\mu$ is no limitation. $\sigma$ is a scale factor so you can set it to $1$ in your initial analysis to see how things work out. A single integration by parts produces a two-term recurrence for $I_m$ in terms of $I_{m-2}$ and $I_{m-1}.$
– whuber
Commented Mar 17 at 15:31
• Yes, I ended up doing the integration by parts! Thank you! Commented Apr 3 at 6:35

Let $$z = (x-\mu)/\sigma$$ and write

$$I_m = \int_0^\infty x^m\, \exp\left(-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2\right)\,\mathrm d x = \int_{-\mu/\sigma}^\infty \left(\sigma z + \mu\right)^m e^{-z^2/2}\,\mathrm dz.$$

Fix $$\mu$$ and $$\sigma.$$ For $$m \ge 1$$ and $$0\le d\lt 1$$ define

$$J_d(m) = \int_{-\mu/\sigma}^\infty (\sigma z + \mu)^{m-1+d}e^{-z^2/2}\,\mathrm dz.$$

An integration by parts gives

\begin{aligned} J_d(m) &= \frac{(\sigma z+\mu)^{m+d}}{\sigma(m+d)}e^{-z^2/2}\bigg|_{-\mu/\sigma}^\infty + \int_{-\mu/\sigma}^\infty \frac{(\sigma z+\mu)^{m+d}}{\sigma(m+d)}\,ze^{-z^2/2}\,\mathrm dz\\ &= \frac{1}{\sigma(m+d)}\int_{-\mu/\sigma}^\infty \frac{1}{\sigma}\left(\sigma z+\mu - \mu\right)\left(\sigma z+\mu\right)^{m+d}e^{-z^2/2}\,\mathrm dz\\ &= \frac{1}{\sigma^2(m+d)}\left(J_d(m+2) - \mu J_d(m+1)\right), \end{aligned}

which upon algebraic rearrangement yields a nice two-term recursion,

$$J_d(m+2) = \sigma^2(m+d)J_d(m) + \mu J_d(m+1).$$

When (as in the question) $$d=0,$$ we may directly compute

$$I_0 = J_0(1) = \sqrt{2\pi}\,\Phi(\mu/\sigma); \quad I_1 = J_0(2) = \sqrt{2\pi}\left(\sigma e^{(\mu/\sigma)^2} + \mu \Phi(\mu/\sigma)\right).$$

to start the recursion. Generally, we could compute $$J_d(1)$$ and $$J_d(2)$$ numerically and then use the recursion to find $$J_d(m)$$ with $$O(m)$$ additional effort.

For concreteness, here is R code for numerical evaluation of $$I_m$$ (for any non-negative real $$m,$$ not just natural numbers) as function f and to implement the recursion (function g). The optional arguments ... are passed to integrate to control its precision and accuracy.

f <- Vectorize(function(m, mu = 0, sigma = 1, ...) {
integrate(\(x) exp(m * log(x) - ((x - mu) / sigma)^2 / 2), 0, Inf, ...)$value }, "m") g <- Vectorize(function(m, mu = 0, sigma = 1, ...) { m.0 <- floor(m) d <- m - floor(m) y <- f(c(d, d + 1), mu, sigma, ...) if (m.0 >= 2) { for (j in seq(1, m.0 - 1)) y[j + 2] = sigma^2 * (j + d) * y[j] + mu * y[j + 1] } head(y, m.0 + 1) }, "m")  As a check and to illustrate, I used this recursion to evaluate $$I_m$$ over a grid for $$\mu\in[-2,3],$$ $$\sigma\in[0.5,2],$$ and $$m=1,2,3,4.$$ The values of f and g are identical to within floating point error. Here are plots of g. • Thank you for your answer! I was just wondering why you would need to introduce$J_d$? I ended up doing an integration by parts to find a recursion on$I_m$directly. Commented Apr 3 at 6:39 • Introducing$d$enabled me to obtain results for all positive values$m$rather than just integral values while, at the same time, showing the recursion still holds. Since it didn't take any more work than otherwise -- only a judicious choice of notation -- I felt having to write a few more "$d\$'s" was worth the cost ;-).
– whuber
Commented Apr 3 at 13:47
• Great, makes sense! Thank you for your input! Commented Apr 5 at 10:17