# How to calculate Estimated Arithmetic Mean for a lognormal distribution

I have been tasked to program functions from some Excel sheet into a asp.net app so it can be shared to my colleagues via a web interface. However, I am stuck on one thing.

I have a set of variables $X = \{0.043, 0.236, 0.057, 0.016\}$ and I have been able to calculate all the descriptive statistics and plot a lognormal graph, etc., but I need to calculate the estimated arithmetic mean which is giving me a headache as I can't find much information on this.

The Excel sheet says the estimated arithmetic mean should be $0.085$. I would be grateful if anybody could sketch me the algorithm needed to get this value as I can't find out how Excel is doing it.

I know it is fit to be normally and lognormally distributed and here are some of the things I have calculated for this series that may help:

Mean: 0.088
Median: 0.05
Standard deviation: 0.1
Geometric mean: 0.0552
Geometric standard deviation: 3.04
W-test of log transformed data: 0.970
W-test of data: 0.786
mean of log(samples): -2.897
standard deviation of log(samples): 1.1117

Edit: Thank you all for your help, it has helped me a lot. I have translated whuber's algorithm in C# code in case anyone needs it:

public double Finney(int m, double z)
{
int i = 0;
int iMax = 0;
double a = 1;
double g = a;
double x = 0;
const double aTol = 0.0000000001;
const int itermax = 1000;

if (m <= -1)
{  return 0;  }

x = (z * m * m) / (m + 1);

if (Math.Abs(x) < aTol)
{  return 1;  }

iMax = Int32.Parse(Math.Round((Math.Abs(Math.Round(z) + 1) + 20)).ToString());

if (iMax > itermax)
{  return 0;  }

for (i = 1; i <= iMax; i++ )
{
if(Math.Abs(a) <= (aTol * Math.Abs(g)))
{  break;  }
a = (a * x) / (m + (2 * (i - 1))) / i;
g = g + a;
}

return g;

}

• Hi and welcome to the site. I really don't know what you mean by "estimated arithmetic mean". The arithmetic mean of the given numbers is $0.088$, as you say. What exactly do you want to calculate? – COOLSerdash Jul 3 '14 at 19:06
• Hey COOLSerdash, Yeah i am confused as well since the description for the estimated arithmetic mean says "Estimated by the Minimum Variance Unbiased Estimate (MVUE), usually more accurate than the simple arithmetic mean of the data" so im not sure what that is, i looked up MVUE but i cannot find anything related to the arithmetic mean – i_shoot_photos Jul 3 '14 at 19:09
• The calculation of the MVUE for the mean of a lognormal distribution is rather complicated and involves infinite series. But check out @whuber's Excel sheet that does the calculations for you. – COOLSerdash Jul 3 '14 at 20:13

## 2 Answers

For positive data $x_1, x_2, \ldots, x_n$ let $y_i = \log(x_i)$ be their natural logarithms. Set

$$\bar{y} = \frac{1}{n}(y_1+y_2+\cdots + y_n)$$

and

$$s^2 = \frac{1}{n-1}\left((y_1 - \bar{y})^2 + \cdots + (y_n - \bar{y})^2\right);$$

these are the mean log and variance of the logs, respectively. The UMVUE for the arithmetic mean when the $x_i$ are assumed to be independent and identically distributed with a common lognormal distribution is given by

$$m(x) = \exp(\bar{y}) g_n\left(\frac{s^2}{2}\right)$$

where $g_n$ is Finney's function

$$g_n(t) = 1 + \frac{(n-1)t}{n} + \frac{(n-1)^3t^2}{2!n^2(n+1)} + \frac{(n-1)^5t^3}{3!n^3(n+1)(n+3)}+\frac{(n-1)^7t^4}{4!n^4(n+1)(n+3)(n+5)} + \cdots.$$

For the data in the question, $s^2 = 1.23594$, $g_4(s^2/2) = 1.532355$, and the UMVUE is $m(x) = 0.084519.$

Because this might take a while to converge when $s^2/2 \gg 1$, it is best implemented as an Excel macro. Such power series are straightforward to program efficiently: just maintain a version of the current term and at each step update it to the next term and add that to a cumulative sum. The term values will typically rise and then fall again; stop when they have fallen below a small positive threshold. (For less floating point error, first compute all such terms and then sum them from smallest to largest in absolute value.)

My version of this macro (in very plain vanilla VBA) follows.

'
' Finney's G (Psi) function as in Millard & Neerchal, formula 5.57
' or equivalently in Gilbert, formula 13.4 (m here = n-1 there).
'
' Typically, m is a positive integer.  Z can be positive or negative.
'
' Programmed by WAH @ QD 5 March 2001
'
' This algorithm will be less accurate for large m*z.  It could be replaced by
' one that separately computes the descending half of the terms,
' iterating backward over i.
'
' It can be badly inaccurate for very negative m*z.
'
' This function returns 0 (an impossible value) upon encountering
' an input error.
'
Public Function Finney(m As Integer, z As Double) As Double
Dim i As Integer    ' Index variable
Dim g As Double     ' Result
Dim x As Double     ' z * m * m / (m+1)
Dim a As Double     ' Power series coefficient
Dim iMax As Integer                 ' Maximum iteration count
Const aTol As Double = 0.0000000001 ' Convergence threshold
Const iterMax As Integer = 1000     ' Limits execution time

If (m <= -1) Then
' issue an error
Finney = 0#
End If

x = z * m * m / (m + 1)

If (Abs(x) < aTol) Then
Finney = 1#     ' This is the correct answer.
Exit Function
End If

iMax = Abs(Int(z) + 1) + 20
If (iMax > iterMax) Then
' issue an error
Finney = 0#
Exit Function
End If
'
' Initialize
'
a = 1#
g = a                       ' Lead terms

For i = 1 To iMax
'
' Test for convergence
'
If (Abs(a) <= aTol * Abs(g)) Then
Exit For
End If
'
' Compute the next term
'
a = a * x / (m + 2 * (i - 1)) / i
'
' Accumulate terms
'
g = g + a
Next

Finney = g
End Function


### References

Gilbert, Richard O. Statistical Methods for Environmental Pollution Monitoring. Van Nostrand Reinhold Company, 1987.

Millard, Steven P. and Nagaraj K. Neerchal, Environmental Statistics with S-Plus. CRC Press, 2001.

### Appendix

For those using a vectorized implementation it pays to precompute the coefficients of $g_n$ in advance for a given value of $n$. This can also be exploited to determine in advance how many coefficients will be needed, thereby avoiding almost all the comparison operations. Here, as an example, is an R implementation. (It uses the equivalent Gamma-function formula of http://www.unc.edu/~haipeng/publication/lnmean.pdf after correcting a typographical error there: the power series argument should be $(n-1)^2t/(2n)$ rather than $(n-1)t/(2n)$ as written.)

finney <- function(t, n, eps=1.0e-20) {
u <- t * (n-1)^2 / (2*n)
tau <- max(u)
i.max <- ceiling(max(1, -log(eps), 1 + log(tau)/2))
a=lgamma((n-1)/2) - (lgamma(1:i.max+1) + lgamma((n-1)/2 + 1:i.max))
b <- exp(a[a + log(tau) * 1:i.max > log(eps)]) # Retain only terms larger than eps
x <- outer(u, 1:length(b), function(z,i) z^i)  # Compute powers of u
return(x %*% b + 1)                            # Sum the power series
}


For example, finney(1.2359357/2, 4) produces the value $1.532355$. This implementation can compute a million values per second for $n=3$ and about $400,000$ values per second for $n=300$. As another example of its use, here is a plot of $g_4, g_8, g_{16}, g_{32}$. (The higher graphs correspond to larger values of $n$.)

par(mfrow=c(1,1))
curve(finney(x/2, 32), 0, 2, lwd=2, main="Finney g(t/2)", xlab="t", ylab="")
curve(finney(x/2, 16), add=TRUE, lwd=2, col="#2040c0")
curve(finney(x/2, 8), add=TRUE, lwd=2, col="#c02040")
curve(finney(x/2, 4), add=TRUE, lwd=2, col="#40c020") • (+1) Great answer. I don't know if this is relevant, but the UMVUE for the arithmetic can also be calculated using the $_0F_1$ Hypergeometric function. Namely $m(x)=\exp{(\bar{y})}_0F_1\left(;\frac{(n-1)}{2};\frac{(n - 1)^{2}s_{y}^{2}}{4n}\right)$. I mention this because the Hypergeometric function may be implemented already (e.g. in R in the hypergeo package) eliminating the need to program Finney's G yourself. – COOLSerdash Jul 3 '14 at 21:35
• @COOL Excellent point! Thank you for sharing that. – whuber Jul 3 '14 at 21:59
• Wow guys thank you so much! This was very useful, and i managed to translate it into c# in case anyone would need the c# version of whuber's algorithm, i will paste it in the original question. – i_shoot_photos Jul 4 '14 at 14:27
• How would one go about using the hypergeometric function to calculate the UMVUE of the variance? – Thomas Jan 12 '15 at 15:43
• @Thomas It's given by $$\exp(2\bar{y}) \left(g_n\left(\frac{2s^2}{n}\right) - g_n \left(\frac{(n-2)s^2}{n(n-1)}\right) \right).$$ – whuber Jan 12 '15 at 15:49

@whuber gave already a complete answer. For convenience, I want to share an implementation of whuber's algorithm in R along with two other solutions using pre-existing packages.

## Using whuber's algorithm

#-----------------------------------------------------------------------------
# The data
#-----------------------------------------------------------------------------

x <- c(0.043, 0.236, 0.057, 0.016)
n <- length(x)
logx <- log(x)

log.mean <- mean(logx)
log.sd <- sd(logx)

#-----------------------------------------------------------------------------
# R-translation of whuber's algorithm "Finney"
#-----------------------------------------------------------------------------

Finney <- function(m, z, maxiter = 1000, aTol = 1e-10){

aTol <- aTol

iterMax <- maxiter

if (m <= -1) {
stop("Finney = 0")
}

x <- z*m*m/(m + 1)

if (abs(x) < aTol) {
return(Finney = 1L)
}

iMax <- abs(trunc(z) + 1) + 20

if (iMax > iterMax) {
stop("iMax > iterMax")
}

a <- 1L
g <- a

for  (i in seq(iMax)) {
if (abs(a) <= aTol*abs(g)) {
break()
}
a <- a*x/(m + 2*(i - 1))/i
g = g + a
}
return(g)
}

# Sanity check

Finney(n-1, log.sd^2/2)
 1.532355

exp(log.mean)*Finney(n-1, log.sd^2/2)
 0.08451876


## Using the hypergeo package

Seems correct. Now the solution using the R package hypergeo. The UMVUE for the arithmetic mean can also be calculated using the $_0F_{1}$ Hypergeometric function in the following way: $$m(x) = \exp{(\bar{y})}_0F_{1}\left(;\frac{(n-1)}{2};\frac{(n-1)^{2}s_{y}^{2}}{4n}\right)$$

#-----------------------------------------------------------------------------
# Using the package "hypergeo"
#-----------------------------------------------------------------------------

require(hypergeo)

genhypergeo(NULL, (n-1)/2, ((n - 1)^2*log.sd^2)/(4*n))
 1.532355

exp(log.mean)*genhypergeo(NULL, (n-1)/2, ((n - 1)^2*log.sd^2)/(4*n))
 0.08451876


## Using the EnvStats package

The package EnvStat has a function elnormAlt that estimates the mean (optionally with a confidence interval) and the coefficient of variation of a lognormal distribution using several methods. Choose the option method = "mvue" to reproduce the results shown above:

#-----------------------------------------------------------------------------
# Using the package "EnvStats"
#-----------------------------------------------------------------------------

require(EnvStats)

elnormAlt(x, method = "mvue", ci = FALSE)

Results of Distribution Parameter Estimation
--------------------------------------------

Assumed Distribution:            Lognormal

Estimated Parameter(s):          mean = 0.08451876
cv   = 1.02389278

Estimation Method:               mvue

Data:                            x

Sample Size:                     4


## Timing the three implementations

Finally, here is a comparison of how long it takes to apply the three methods to 1,000 samples of size $n=5,10,15,...,1000$, using @whuber's method as the baseline. The functions from the EnvStats and hypergeo packages presumably have more error handling and more options, which at least partially can explain why they take so much longer. The R code used for the comparison follows below:

nvec<-seq(10,1000,10)
B<-1000
reftime<-time1<-time2<-time3<-rep(NA,length(nvec))

# Compile the COOLSerdash-Whuber function:
require(compiler)
Finney<-cmpfun(Finney)

for(i in 1:length(nvec))
{
n<-nvec[i]
cat(n,"\n")

## Just generate some LNorm data:
start.time <- Sys.time()
for(j in 1:B) {x<-rlnorm(n)}
end.time <- Sys.time()
reftime[i]<-end.time - start.time

## Whuber's method:
start.time <- Sys.time()
for(j in 1:B) {x<-rlnorm(n); exp(log.mean)*Finney(n-1, log.sd^2/2)
}
end.time <- Sys.time()
time1[i]<-end.time - start.time

## Hypergeo:
start.time <- Sys.time()
for(j in 1:B) {x<-rlnorm(n); exp(log.mean)*genhypergeo(NULL, (n-1)/2, ((n - 1)^2*log.sd^2)/(4*n))
}
end.time <- Sys.time()
time2[i]<-end.time - start.time

## EnvStats:
start.time <- Sys.time()
for(j in 1:B) {x<-rlnorm(n); elnormAlt(x, method = "mvue", ci = FALSE)  }
end.time <- Sys.time()
time3[i]<-end.time - start.time

}

##
time1<-time1-reftime
time2<-time2-reftime
time3<-time3-reftime

## Save the results:
plot(nvec,time1,type="l",lwd=3,ylim=c(0,max(time3/time1)),ylab="Relative execution time",xlab="Sample size n",cex.lab=1.5,cex.axis=1.5,cex.main=1.5,main="Relative execution time")
lines(nvec,time2/time1,type="l",lwd=3,col=2)
lines(nvec,time3/time1,type="l",lwd=3,col=4)
legend(600,19,c("@whuber","hypergeo","EnvStats"),col=c(1,2,4),lwd=2,cex=1.5)