I am trying to solve the integral,
\begin{equation} \int_{0}^b\Phi_{LN}(x, \mu,\sigma)dx \end{equation}
where $\Phi_{LN}(x, \mu,\sigma)$ is the lognormal cdf, evaluated at $x$. i.e., $$\int^x_{0}\frac{1}{\sqrt{2\pi\sigma}y}\exp\left({-\frac{1}{2}\left(\frac{\log y-\mu}{\sigma}\right)^2}\right)dy$$.
Being able to efficiently compute the outer integral (top equation) would be great. This would ideally involve expressing it in a form that most software could efficiently solve/approximate.
I thought I could have solved this problem by solving the integral to a standard normal cdf, and simply exchanging $x$ with $\frac{log(x)-\mu}{\sigma}$, however this does not seemed to have worked.
In case someone wants to use the integral of the standard normal cdf to answer this question it is shown following,
$$ \int\Phi_{N(0,1)}(x) = x\Phi_{N(0,1)}(x) - \int x\Phi^{'}_{N(0,1)}(x)dx $$
and here is my failed attempt to convert it to work as the lognormal equivalent,
> # define function to integrate N(0,1) CDF
> int_LN_CDF = function(x, mu, sigma){
+ # transform to lognormal RV
+ y = (log(x) - mu)/sigma
+ # integral of standard normal cdf
+ y*pnorm(y) + 1/sqrt(2*pi)*exp(-(y^2)/2)
+ }
> # set arbitary values of mu and sigma
> mu = log(2); sigma = 0.4
> # test function against numerical approach
> int_LN_CDF(2, mu, sigma) - int_LN_CDF(1, mu, sigma)
[1] 0.3820694
> integrate(function(y) pnorm(y, mu, sigma), lower = 1, upper = 2)
[1] 0.9491263 with absolute error < 1.1e-14
Any help or guidance would be appreciated.
Thanks
You can use a similar approach as for the standard normal distribution used here solving/approximating integral of standard normal cdf. I'll denote by $\Phi$ and $\varphi$ the distribution and density function of the standard normal distribution - the general case should be feasible from there. The distribution function of the log-normal distribution is then $F(x) = \Phi(\log(x))$. Starting with \begin{align} \frac{d}{dx} xF(x) = F(x) + \varphi(\log(x)) \end{align} we need to find a function $L(x)$ such that $\frac{d}{dx} L(x) = -\varphi(\log(x))$. Then we can finally define $I(x) = xF(x) + L(x)$ with $\frac{d}{dx} I(x) = F(x)$ which is what you are looking for. Finding $L$ is not straight forward - so I asked Wolfram Alpha for help. With $L(x) = \frac{1}{2 e^{-0.5}} \mbox{erf} ((1-\log(x)) / \sqrt{2})$, where $\mbox{erf}$ is the error function \begin{align} \mbox{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} dt \end{align} you should get your result.
Edit: to make my suggestion in the comment concerning the general case more explicit I suggest to look at \begin{align} L(x) = \frac{1}{2}e^{(\mu+\sigma^2/2)}\mbox{erf}\left(\left(\sigma - \frac{\log(x)-\mu}{\sigma}\right)/\sqrt{2}\right). \end{align} With $\frac{d}{dx}\mbox{erf}(x) = \frac{2}{\sqrt{\pi}}e^{-x²}$ we then get \begin{align} \frac{d}{dx} L(x) &= -\frac{1}{2}e^{(\mu+\sigma^2/2)}\frac{2}{\sqrt{\pi}}e^{-\frac{1}{2}\left(\sigma - \frac{\log(x)-\mu}{\sigma}\right)^2}\frac{1}{x\sigma\sqrt{2}}\\ &= -\frac{1}{x\sigma\sqrt{2\pi}}e^{(\mu+\sigma^2/2)}e^{-\frac{1}{2}\left(\sigma^2 + \left(\frac{\log(x)-\mu}{\sigma}\right)^2 - 2(\log(x)-\mu)\right)}\\ &= -\frac{1}{x\sigma\sqrt{2\pi}}e^{(\mu+\sigma^2/2)}e^{-\frac{1}{2}\sigma^2}e^{-\frac{1}{2}\left(\frac{\log(x)-\mu}{\sigma}\right)^2}e^{\log(x)}e^{-\mu}\\ &= -\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{\log(x)-\mu}{\sigma}\right)^2}\\ &= -\frac{1}{\sigma}\varphi\left(\frac{\log(x)-\mu}{\sigma}\right). \end{align} Turning back to the general problem with $F_{\mu,\sigma}(x) = \Phi\left(\frac{\log(x)-\mu}{\sigma}\right)$ we now go on defining $I(x) = xF_{\mu,\sigma}(x) + L(x)$ to get \begin{align} \frac{d}{dx}I(x) &= F_{\mu,\sigma}(x) + xf_{\mu,\sigma}(x) - \frac{1}{\sigma }\varphi\left(\frac{\log(x)-\mu}{\sigma}\right)\\ &= F_{\mu,\sigma}(x), \end{align} i.e. $I(x)$ is the anti-derivative of $F_{\mu,\sigma}(x)$.