Trying to solve a homework problem:
I have two functions to get the geometric mean from 1000 observations from the exponential distribution with a rate of .01. The following keeps returning Inf.
gmean <- function(n)
{
prod(n)^(1/length(n))
}
x<-rexp(1000,1/100)
gmean(x)
but this does not
gmean1 <- function(n)
{
exp(mean(log(n)))
}
x<-rexp(1000,1/100)
gmean1(x)
Why is this? I think it's something to do with the prod function but I'm not sure.
The problem is that when you do prod(n) in your function, it calculates the result of this call before raising it to the power of (1/length(n)). Since the mean of x is about 100, you can expect this call to return a value with a similar order of magnitude to 100^1000, which is much higher than the maximum number that R will return (R will call anything above around 10^308 Inf).
Any mathematical operation you attempt on Inf will also return Inf, so your naive implementation will not work if x is greater than about 154:
100^154
#> [1] 1e+308
100^155
#> [1] Inf
In actuality, because the majority of numbers are less than 100 in your sample, you might get to an x length of about 180 before you started generating Inf
In any case, it would be safer to stick to
gmean <- function(n) exp(sum(log(n))/length(n))
