Calculate mean of density function (example using lognormal distribution)

In a lognormal distribution, the mean is equal to $$\exp(\mu + \frac{\sigma^2}{2})$$.

I tried to separately calculate this using the definition $$E[X] = \int{xf(x)dx}$$, where I have 200 $$x$$ values evenly log-spaced; for reproducibility, here a few copy/paste ways to get such an array that captures the tails when $$\mu=-0.5$$ and $$\sigma=0.6$$:

• in R: x <- exp(seq(log(0.001), log(130), length.out = 200))
• in Python: x = np.logspace(start=np.log(0.001), stop=np.log(130), num=200, base=np.e)
• in Matlab: x = exp((linspace(log(.001), log(130), 200)))';

In Excel, I have the following steps based off the x-array defined above:

1. Use the NORM.DIST function with x=ln(x), mean=-0.5, standard_dev=0.6, cumulative=FALSE to get the density function $$f(x)$$

2. Since the $$f(x)dx$$ term in the $$E[X]$$ integral is just the area, I calculate the area for each step using the trapezoidal rule: $$area = \frac{f(x_i) + f(x_{i+1})}{2} \times (x_{i} - x_{i+1})$$

3. Calculate weighted average using SUMPRODUCT function with array1=x, array2=area.

#3 results in 0.735142709; but directly calculating it from the definition of mean for lognormal, the answer should be 0.726149037. I think this x-array is sufficiently dense for the answers to be closer. If fact, if you just treat the density column as a PMF instead of a PDF, you get the exact correct answer (i.e., =(SUMPRODUCT=array1=x, array2=area)/sum(array=density)).

What am I missing that is causing #3 to not give the exact answer, when I essentially have PMF here? (Note that when the area is computed in semi-log space--i.e., $$area = \frac{f(x_i) + f(x_{i+1})}{2} \times (ln(x_{i}) - ln(x_{i+1}))$$-- it sums to unity; using this area in the integral doesn't work either).

The overall motivation for this question is I am trying to calculate the mean for a distribution that doesn't have a closed-form solution. I tried this approach here as a proof-of-concept, expecting it to work, so I am trying to figure out what I am misunderstanding. Thank you for any help.

Copy/paste for Excel:

If you have the x-array in column A, then:

• ln(x) in column B =LN(A2)
• density in column C =NORM.DIST(B2,-0.5,0.6,FALSE)
• area[arith] in column D =(C2+C3)/2*(A3-A2)
• application of integral is =SUMPRODUCT(A2:A201,D2:D201)
• area[ln] in column E is =(C2+C3)/2*(B3-B2)

I'm not exactly sure what is your code, but your description does not sound right.

Trapezoidal rule tells us that we can approximate the integral with a sum

$$\int f(x) \, dx \approx \sum_i \tfrac{f(x_{i+1}) + f(x_i)}{2} (x_{i+1} - x_i)$$

What you seem to be saying is that you calculated

$$area_i = \tfrac{f(x_{i+1}) + f(x_i)}{2} (x_{i+1} - x_i)$$

and use it to get

$$\sum_i area_i \times x_i$$

The problem is that this is that

$$\frac{f(x_{i+1}) + f(x_i)}{2} (x_{i+1} - x_i) x_i \\ \ne \frac{x_{i+1} f(x_{i+1}) + x_i f(x_i)}{2} (x_{i+1} - x_i)$$

as you are trying to approximate

$$\int x \, f(x) \, dx$$

Also keep in mind that this is just an approximation, and a rather rudimentary one. With only 200 points you should not expect it to be very precise, even if implemented correctly. If you want to do numerical integration, there are many better solutions, available out of the box, like the integrate function in R, the scipy.integrate subpackage in Scipy, etc. Excel is usually a rather poor solution for such problems.

• Why does =(SUMPRODUCT=array1=x, array2=area)/sum(array=density) return exactly the correct answer?
– a11
Jul 10, 2023 at 22:11
• @a11 That's a different question and a longer answer, but your solution is incorrect, and because of that slightly off (it uses $x_i$ in a place of $x_{i+1}$ if you look carefully). The alternative approach accidentally is closer to a different approximation of the integral.
– Tim
Jul 10, 2023 at 22:20
• "area[arith] in column D =(C2+C3)/2*(A3-A2)" and "application of integral is =SUMPRODUCT(A2:A201,D2:D201)" column A is x, so it is using x, not x(i+1)?
– a11
Jul 10, 2023 at 22:23
• regardless, your comment about using something other than Excel was a good point; the trapz function in R pracma produces the expected answer: trapz(x,density)
– a11
Jul 10, 2023 at 22:30
• @a11 I'd think of it rather as of histogram approximation of pdf than as of pmf. You are taking a weighted average there, where weights are probability densities, but for a weighted average, the weights are normalized so they sum to 1.
– Tim
Jul 11, 2023 at 7:32
1. In principle, the numeric integration looks almost correct. To get better results, you should however take the sum over midpoint.x times area, where the $$i$$-th midpoint is $$(x_i+ x_{i+1})/2$$. What you then calculate, is $$\sum_{i=1}^n \frac{x_{i+1}+x_{i}}{2}\frac{f(x_i)+f(x_{i+1})}{2}{(x_{i+1}-{x_i})}$$

2. The density of the lognormal distribution is not calculated correctly, it should read $$f(x)=\frac{1}{x}\phi_{\mu,\sigma}(\ln(x))$$, where $$\phi_{\mu,\sigma}$$ is the normal density.

in Excel:

• replace formula for column C by =NORM.DIST(B2, -0.5, 0.6, FALSE)/A2
• add column F = (A2 + A3)/2 then the integral is =SUMPRODUCT(D2:D201,F2:D201). This will give you the result 0.72783.

The way to approximate the integral by multiplying the trapezoidal areas with the midpoints is a kind of hybrid between rectangular and trapezoidal rule.

For a proper implementation of the trapezoidal rule on the integrand $$x f(x)$$ in excel, you can use that $$xf(x) = \phi_{\mu,\sigma}(\ln(x))$$, generate H=NORM.DIST(B2;-0,5; 0,6;FALSE), and I =(H3+H2)/2*(A3-A2), and calculate the integral as sum(I2:I201). The result is then 0.726569. This trick works only for special cases as the log normal distribution, but not in general.

To implement the trapezoidal rule in the general case, assume you have the correct density $$f(x)$$ in column C and the $$x$$ values in A. Then create D =(C3*A3+C2*A2)/2*(A3-A2) and calculate the integral as =sum(D2:D201).

As @Tim wrote, numerical integration is easy and comfortable in R or python, but maybe you are bound to do this in excel.

Another reason for discrepancy between the theoretical result and numeric integration are the finite limits for the integral. With R's function integrate, you can give infinite upper and lower limits, and specify the accuracy of the numeric integration.

> integrate(function(x) dlnorm(x, -.5, .6) * x,
lower = 0, upper=Inf, rel.tol = 0.000001)
0.726149 with absolute error < 5.6e-08