# Is it possible to analytically integrate $x$ multiplied by the lognormal probability density function?

Firstly, by analytically integrate, I mean, is there an integration rule to solve this as opposed to numerical analyses (such as trapezoidal, Gauss-Legendre or Simpson's rules)?

I have a function $\newcommand{\rd}{\mathrm{d}}f(x) = x g(x; \mu, \sigma)$ where $$g(x; \mu, \sigma) = \frac{1}{\sigma x \sqrt{2\pi}} e^{-\frac{1}{2\sigma^2}(\log(x) - \mu)^2}$$ is the probability density function of a lognormal distribution with parameters $\mu$ and $\sigma$. Below, I'll abbreviate the notation to $g(x)$ and use $G(x)$ for the cumulative distribution function.

I need to calculate the integral $$\int_{a}^{b} f(x) \,\rd x \>.$$

Currently, I'm doing this with numerical integration using the Gauss-Legendre method. Because I need to run this a large number of times, performance is important. Before I look into optimizing the numerical analyses/other pieces, I would like to know if there are any integration rules to solve this.

I tried applying the integration-by-parts rule, and I got to this, where I'm stuck again,

1. $\int u \,\mathrm{d}v = u v - \int v \mathrm{d}u$.

2. $u=x \implies \rd u = \rd x$

3. $\rd v = g(x) \rd x \implies v = G(x)$

4. $u v - \int v \rd x = x G(x) - \int G(x) \rd x$

I'm stuck, as I can't evaluate the $\int G(x) \rd x$.

This is for a software package I'm building.

• @Rosh, by $lognormal$ you mean probability density of log-normal distribution? – mpiktas Apr 13 '11 at 12:45
• This is expressible as a constant times a difference of two normal cdfs. Normal cdfs are efficiently computed using W. Cody's rational Chebyshev approximation. You shouldn't need and, almost undoubtedly shouldn't prefer, numerical-integration alternatives to this. If you need more details, I can post them. – cardinal Apr 13 '11 at 12:50
• @mpiktas, Yes, lognormal is the probability density function and lognormalCDF is the cumulative density function. – Rosh Apr 13 '11 at 13:23
• @Rosh $x$ has a lognormal distribution means that $\log(x)$ is normally distributed. Thus, substitute $x = \exp(y)$ in your original integral. The integrand is an exponential whose argument is a quadratic function of $y$. Completing the square turns it into a multiple of a normal PDF, so your answer is written in terms of the normal CDF and the exponentials of the original endpoints. There are many good approximations to the normal CDF (a multiple of the error function). – whuber Apr 13 '11 at 14:12
• Yes, @whuber and I were describing the same thing. You should get something like $e^{\mu + \frac{1}{2} \sigma^2} (\Phi(\beta) - \Phi(\alpha))$ where $\beta = (\log(b) - (\mu + \sigma^2))/\sigma$ and $\alpha = (\log(a) - (\mu+\sigma^2))/\sigma$ and $\Phi(\cdot)$ denotes the normal cdf. Note that, depending on the values of $a$, $b$, $\mu$ and $\sigma$, there are ways to rewrite this expression to be more numerically stable. – cardinal Apr 13 '11 at 17:19

Short answer: No, it is not possible, at least in terms of elementary functions. However, very good (and reasonably fast!) numerical algorithms exist to calculate such a quantity and they should be preferred over any numerical integration technique in this case.

Quantity of interest in terms of normal cdf

The quantity you are interested in is actually closely related to the conditional mean of a lognormal random variable. That is, if $X$ is distributed as a lognormal with parameters $\mu$ and $\sigma$, then, using your notation, $$\newcommand{\e}{\mathbb{E}}\renewcommand{\Pr}{\mathbb{P}}\newcommand{\rd}{\mathrm{d}} \int_a^b f(x) \rd x = \int_a^b \frac{1}{\sigma \sqrt{2\pi}} e^{-\frac{1}{2\sigma^2}(\log(x) - \mu)^2} \rd x = \Pr(a \leq X \leq b) \e(X \mid a \leq X \leq b) \>.$$

To get an expression for this integral, make the substitution $z = (\log(x) - (\mu + \sigma^2))/\sigma$. This may at first appear a bit unmotivated. But, note that using this substitution, $x = e^{\mu + \sigma^2} e^{\sigma z}$ and by a simply change of variables, we get $$\int_a^b f(x) \rd x = e^{\mu + \frac{1}{2}\sigma^2} \int_{\alpha}^{\beta} \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2} z^2} \rd z \> ,$$ where $\alpha = (\log(a) - (\mu + \sigma^2))/\sigma$ and $\beta = (\log(b) - (\mu + \sigma^2))/\sigma$.

Hence, $$\int_a^b f(x) \rd x = e^{\mu + \frac{1}{2}\sigma^2} \big( \Phi(\beta) - \Phi(\alpha) \big) \>,$$ where $\Phi(x) = \int_{-\infty}^x \frac{1}{\sqrt{2\pi}} e^{-z^2/2} \rd z$ is the standard normal cumulative distribution function.

Numerical approximation

It is often stated that no known closed form expression for $\Phi(x)$ exists. However, a theorem of Liouville from the early 1800's asserts something stronger: There is no closed form expression for this function. (For the proof in this particular case, see Brian Conrad's writeup.)

Thus, we are left to use a numerical algorithm to approximate the desired quantity. This can be done to within IEEE double-precision floating point via an algorithm of W. J. Cody's. It is the standard algorithm for this problem, and utilizing rational expressions of a fairly low order, it's pretty efficient, too.

Here is a reference that discusses the approximation:

W. J. Cody, Rational Chebyshev Approximations for the Error Function, Math. Comp., 1969, pp. 631--637.

It is also the implementation used in both MATLAB and $R$, among others, in case those make it easier to obtain example code.

Here is a related question, in case you're interested.