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;

  • 2
    $\begingroup$ 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? $\endgroup$ Jul 3, 2014 at 19:06
  • 2
    $\begingroup$ 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 $\endgroup$ Jul 3, 2014 at 19:09
  • 3
    $\begingroup$ 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. $\endgroup$ Jul 3, 2014 at 20:13

2 Answers 2


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)$$


$$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

    Finney = g
End Function


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.


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$.)

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")


  • 4
    $\begingroup$ (+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. $\endgroup$ Jul 3, 2014 at 21:35
  • 2
    $\begingroup$ @COOL Excellent point! Thank you for sharing that. $\endgroup$
    – whuber
    Jul 3, 2014 at 21:59
  • $\begingroup$ 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. $\endgroup$ Jul 4, 2014 at 14:27
  • $\begingroup$ How would one go about using the hypergeometric function to calculate the UMVUE of the variance? $\endgroup$
    – Thomas
    Jan 12, 2015 at 15:43
  • $\begingroup$ @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).$$ $\endgroup$
    – whuber
    Jan 12, 2015 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)) {
    a <- a*x/(m + 2*(i - 1))/i
    g = g + a

# Sanity check

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

exp(log.mean)*Finney(n-1, log.sd^2/2)
[1] 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"


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

exp(log.mean)*genhypergeo(NULL, (n-1)/2, ((n - 1)^2*log.sd^2)/(4*n))
[1] 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"


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.

enter image description here

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:


# Compile the COOLSerdash-Whuber function:

for(i in 1:length(nvec))

## 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



## 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")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.