Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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;

    }
share|improve this question
2  
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 at 19:06
2  
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 at 19:09
3  
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 at 20:13

2 Answers 2

up vote 5 down vote accepted

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

Graphs

share|improve this answer
2  
(+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 at 21:35
1  
@COOL Excellent point! Thank you for sharing that. –  whuber Jul 3 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 at 14:27

@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] 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"
#-----------------------------------------------------------------------------

require(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"
#-----------------------------------------------------------------------------

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.

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:

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)
share|improve this answer
    
Very nice (+1), @COOL! Just for fun I ran some simulations comparing the execution times of the three implementations. Do you mind if I add a figure describing the results to your post? –  MånsT Jul 4 at 7:26
    
@MånsT That's cool! No, I welcome additions to my posts, feel free to add the figure. –  COOLSerdash Jul 4 at 7:31
    
@MånsT Very nice! Thank you for your effort! –  COOLSerdash Jul 4 at 7:49
    
Thank you very much MansT and COOLSerdash :) –  i_shoot_photos Jul 4 at 14:31

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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