# Approximating lognormal sum pdf (in R)

I have an application for which I need an approximation to the lognormal sum pdf for use as part of a likelihood function. The lognormal sum distribution has no closed form, and there are a bunch of papers in signal processing journals about different approximations. I have been using one of the simplest approximations (Fenton 1960), which involves replacing a sum of lognormals with a single lognormal with matching first and second moments. This is pretty straightforward to code, but judging by the literature on the subject that has been written in the last 50 years, this may not be the best approximation for all applications. I have no intuition for how to identify which approximations will lead to the best MLE estimates.

Does anyone know if (A) There is different approximation I should be using for a maximum likelihood application? (B) There is existing R code for any of the more computationally intensive approximations?

Update: For some background on the problem, see this review

• Can you clarify just a touch? Is what you refer to as the "lognormal sum pdf" the density function of $Y = X_1 + \cdots + X_n$ where $X_n$ are iid lognormal with parameters $\mu$ and $\sigma^2$? – cardinal Jun 4 '11 at 15:10
• Yes, the pdf for the sum of N iid lognormal variates. – Ben Lauderdale Jun 4 '11 at 15:13
• How large is $n$ in your application? – cardinal Jun 4 '11 at 15:14
• I am most interested in the cases where N is small, < 10 or so. However, it would be very helpful if I could at least manage N up to 100 or so. – Ben Lauderdale Jun 4 '11 at 15:16
• Moment matching a lognormal to this sounds on the surface like a strange idea. This is because the lognormal is not characterized by its moments. I will look here, but perhaps there is a way to turn the problem around a little. Let $f_0(x)$ be a "standard" ($\mu = 0$, $\sigma = 1$) lognormal density. For $b \in (-1,1)$, define $f_b(x) = f_0(x) (1 + b \sin(2\pi\log x))$. Then $f_b$ is a pdf and $f_0$ and $f_b$ have the same moments for every such $b$. – cardinal Jun 4 '11 at 15:22

To obtain a numerical version of distribution function for moderate $N$ (say a dozen of r.vs or less), a simple approach is to compute the Discrete Fourier Transform (DFT) of each LN density, form the product and then use inverse DFT. The same grid must be used for all densities and it must be designed with some care. The computation can be done quite easily in a R function. However, do not expect to reach the remarkable precision of the classical distributions functions in R.