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$– NebulousRevealJan 8, 2011 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 TarczynskiJan 8, 2011 at 13:39
4 Answers
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$.
(Added)
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$ Jan 13, 2011 at 21:25
-
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$ Apr 7, 2011 at 5:48
-
1$\begingroup$ The link to
sciencedirect.com
is broken. I'm also unable to find any snapshot saved on the Wayback Machine. $\endgroup$ May 8 at 17:49 -
1$\begingroup$ @TheAmplitwist I'm sorry about that. I didn't keep records of the results of that search. A new search for "polynomial approximation Bessel function" returns some classics as well as more recent papers such as Li, Li, and Gross on ScienceDirect. (That's on the J function, not I, but the existence of such hits is promising...) $\endgroup$– whuberMay 8 at 18:14
-
1$\begingroup$ @whuber Ah, quite understandable, it was over a decade ago after all. Thanks for taking a look, and for the reference to the paper by Li, Li, and Gross. :) $\endgroup$ May 8 at 18:17
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$ Jan 8, 2011 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$ Jan 9, 2011 at 21:34
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$ Jan 8, 2011 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$ Jan 8, 2011 at 23:31