7
$\begingroup$

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...

Thanks in advance!

$\endgroup$
3
  • $\begingroup$ 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$ $\endgroup$
    – Henry
    Commented Mar 17 at 12:32
  • 1
    $\begingroup$ 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}.$ $\endgroup$
    – whuber
    Commented Mar 17 at 15:31
  • $\begingroup$ Yes, I ended up doing the integration by parts! Thank you! $\endgroup$ Commented Apr 3 at 6:35

1 Answer 1

8
$\begingroup$

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.

enter image description here

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ Commented Apr 3 at 6:39
  • $\begingroup$ 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 ;-). $\endgroup$
    – whuber
    Commented Apr 3 at 13:47
  • $\begingroup$ Great, makes sense! Thank you for your input! $\endgroup$ Commented Apr 5 at 10:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.