Is there any good way to approximate following integral?
$$\int_0^{0.5}\frac{x^2}{\sqrt{2\pi}\sigma}\cdot \exp\left(-\frac{(x^2-\mu)^2}{2\sigma^2}\right)\mathrm dx$$
$\mu$ is between $0$ and $0.25$, the problem is in $\sigma$ which is always positive, but it can be arbitrarily small.
I was trying to expand it using Taylor series, but terms looks more or less this $\pm a_n\cdot\frac{x^{2n+3}}{\sigma^{2n}}$ and that can be arbitrarily large, so the error is significant.

  • $\begingroup$ looks similar to the normal distribution. $\endgroup$ – PrimeNumber Jan 8 '11 at 13:36
  • $\begingroup$ @Trevor: Indeed, it is normal distribution with some 'minor' changes. It is multiplied by $x^2$ and in exponent there is $x^{2}$ instead of $x$. $\endgroup$ – Tomek Tarczynski Jan 8 '11 at 13:39

If you write y=x^2 and pull the constants out you have $$\frac{1}{2\sqrt{2\pi}\sigma}\int_0^{0.25}\sqrt{y}\cdot \exp(-\frac{(y-\mu )^2}{2\sigma ^2})dy$$ If $\sigma$ is very small, the contribution will all come from a small area in $y$ around $\mu$. So you can set $\sqrt{y}=\sqrt{\mu}$ and use your error function tables for a close approximation. A quick search didn't turn up moments of $\sqrt{y}$ against the normal distribution, but maybe they are out there.

  • $\begingroup$ Simple and clever, thanks! I've tested it and it seems that this approximation is close enough. I will wait one more day before accepting Your answer. $\endgroup$ – Tomek Tarczynski Jan 8 '11 at 22:08
  • $\begingroup$ @Tomek: I note that when I put dy=xdx I dropped a factor 2. If you found the approx good enough, probably you found this already. I'll fix. $\endgroup$ – Ross Millikan Jan 9 '11 at 21:34

A standard way to get a good approximation for integrals that "look" Gaussian is to evaluate the Taylor series of the logarithms of their integrands through second order, expanding around the point of maximum value thus (continuing with @Ross Millikan's substitution):

$$\eqalign{ &\log\left(\sqrt{y}\cdot \exp\left(-\frac{(y-\mu )^2}{2\sigma ^2}\right)\right) \cr = &\frac{-\mu ^2-\sigma ^2+\mu \sqrt{\mu ^2+2 \sigma ^2}+2 \sigma ^2 \log\left[\frac{1}{2} \left(\mu +\sqrt{\mu ^2+2 \sigma ^2}\right)\right]}{4 \sigma ^2} \cr + &\left(-\frac{1}{2 \sigma ^2}-\frac{1}{\left(\mu +\sqrt{\mu ^2+2 \sigma ^2}\right)^2}\right) \left(y-\frac{1}{2} \left(\mu +\sqrt{\mu ^2+2 \sigma ^2}\right)\right)^2 \cr + &O\left[y-\frac{1}{2} \left(\mu +\sqrt{\mu ^2+2 \sigma ^2}\right)\right]^3 \cr \equiv &\log(C) - (y - \nu)^2/(2\tau^2)\text{,} }$$

say, with the parameters $C$, $\nu$, and $\tau$ depending on $\mu$ and $\sigma$ as you can see. The resulting integral now is a Gaussian, which can be computed (or approximated or looked up) in the usual ways. The approximation is superb for small $\sigma$ or large $\mu$ and still ok otherwise.

The plot shows the original integrand in red (dashed), this approximation in blue, and the simpler approximation afforded by replacing $\sqrt{y} \to \sqrt{\mu}$ in gold for $\sigma = \mu = 1/20$.

alt text


Mathematica tells us the integral, when taken to $\infty$, can be expressed as a linear combination of modified Bessel Functions $I_\nu$ of orders $\nu = -1/4, 1/4, 3/4, 5/4$ with common argument $\mu^2/(4 \sigma^2)$. From the Taylor expansion we can see that when both $\mu$ and $\sigma$ are small w.r.t. $1/2$--specifically, $(1/4-\mu)/\sigma \gg 3$, the error made by including the entire right tail will be very small. (With a little algebra and some simple estimates we can even get good explicit bounds on the error as a function of $\mu$ and $\sigma$.) There are many ways to compute or approximate Bessel functions, including polynomial approximations. From looking at graphs of the integrand, it appears that the cases where the Bessel function approximation works extremely well more or less complement the cases where the preceding "saddlepoint approximation" works extremely well.

  • $\begingroup$ Great answer, taking logarithm removes the problem with $\sigma$ in the denominator. The second thing I learned from this is that it is better to approximate around real maximum, instead of approximating around approximated maximum (I was calculating this series around $\mu$ to have simpler calculations). $\endgroup$ – Tomek Tarczynski Jan 13 '11 at 21:25
  • $\begingroup$ @Tomek Thanks. There really is no problem with $\sigma$ in the denominator: it's just a scale factor that is easily made to go away with an appropriate substitution. The difficulty lies in $\sqrt{y}$ multiplying a Gaussian involving $y-\mu$ as the variable. You have the correct insight: the approximation usually works best when we expand around the maximum value of the exponent. You can see this is a general approach; it's closely related to saddlepoint approximation and the method of stationary phase. $\endgroup$ – whuber Jan 13 '11 at 21:35
  • $\begingroup$ If we didn't logarithm this function then we would have sigma in the denominator of every derivative and that was the thing that stopped me from using this kind of approximation. $\endgroup$ – Tomek Tarczynski Jan 13 '11 at 21:58
  • $\begingroup$ @Tomek The substitution $y \to z\sigma $ makes $\sigma$ disappear (except in trivial ways as a multiplicative factor and in the upper limit of the integral). This changes the value of $\mu$, but so what... $\endgroup$ – whuber Jan 13 '11 at 22:03
  • 1
    $\begingroup$ I'll just add the note that (often,) quarter-order Bessel functions are more profitably expressed as parabolic cylinder functions (or alternatively, Hermite functions). $\endgroup$ – J. M. is a poor mathematician Apr 7 '11 at 5:48

Not an answer, but might still be helpful. Using the variable substitution that Ross mentions in his answer, we can treat a simpler case $\mu=0$ more easily.

For the following integral, Wolfram Alpha tells us that (I hope I did not make a transcription error here):

$$\int_0^\infty \sqrt{y}e^{-y^2/(2\sigma^2)}dy = \frac{\sigma^{3/2}\Gamma(3/4)}{2^{1/4}},$$

But your problem goes from $0$ to $0.25$, so some approximations are needed.

Link-text: http://www.wolframalpha.com/input/?i=Integrate[Sqrt[x]+Exp[-x^2%2F%282+s^2%29]%2C{x%2C0%2Cinf}]

  • $\begingroup$ Your solution might be useful when $\mu$ is close to $0$, but unfortunately it is not guaranteed. Still +1. $\endgroup$ – Tomek Tarczynski Jan 8 '11 at 22:11

How about some good old-fashioned trapezoid rule?

  • $\begingroup$ I was thinking about it since it is quite easy to predict where the highest value is, it will be somewhere near $\mu$. Due to $3\sigma$ law we know on what interval values are concentrated. We could choose following intervals $[\mu-2\sigma,\mu-\sigma],[\mu-\sigma,\mu],[\mu,\mu+\sigma],[\mu+\sigma,\mu+2\sigma]$, but there is a problem when one of this intervals is outside the boundries. This problem is likely to be solved, nevertheless this approximation would probably be not so 'short'. $\endgroup$ – Tomek Tarczynski Jan 8 '11 at 23:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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