Generating a Random Value Vector from an Exponential Distribution using R

Question

Given a standard PDF of the form $f(x)=ae^{-ax}$ with domain $[0,+\infty)$, its CDF being $F(x)=1-e^{-ax}$, and a mutated CDF that takes $p \in [0,1]$ as a probability and returns the corresponding $x$ value at that percentile: $x=-\frac{ln(1-p)}{a}$ (quantile equation)

I am trying to write an R program that returns a randomized vector of size n of from this PDF distribution. The function below represents the quantile equation from above.

makeExpDist <- function(a, p) {
  return(log(1-p)/(-a))
}

I then use this last function to map a randomized vector (size n) from a uniformly distributed variable $[0,1]$ into the above function, expecting that the returned vector will represent a sample from the original PDF.

giveSample <- function(cdf, a, n){
  random <- runif(n)
  le <- c()
  for(i in random){
    le <- append(le, cdf(a,i))
  }
  return(data.frame(id = 1:n, le_prob = random, le_num = le))
}

I would expect the mean to be $86.6434$ (rounded) since $-\frac{ln(.5)}{.008}=86.6434$ (rounded), which can be supported by $\int_{0}^{86.64}.008e^{-.008x}dx=.5$

To demonstrate:

makeExpDist(.008,.5)

returns

[1] 86.6434

I would also expect that to be the mean of a sufficiently large sample from the PDF. I try this below, but the mean of the resulting vector is always $~125$:

testframe <- giveSample(makeExpDist, .008, 100000) #test
mean(testframe$le_num)
mean(testframe$le_prob)

... which on my last run printed:

[1] 125.7898
[1] 0.5016412

The value of the second mean as expected is around 0$.5$ since that simply is the uniform variable I used.

So my question is why do I keep receiving $~125$ as the mean of randomized values when the true mean of the PDF is $~86$? I have tried many fixes in R and get the same exact mean of $~125$. Also, when I run the same process in excel of mapping a randomized uniform vector onto this distribution, I get the same output of $~125$. Is it wrong to use a random uniform variable $[0,1]$ to randomize a value in the PDF? Given the expected value of runif() is 0.5 why is the expected value of the return vector not equal to makeExpDist(.008,.5)?

Thanks a bunch!

Answer 1

The quantity you calculated (86.6434) is the 50th percentile, ie the median, not the mean. Since the exponential is a skewed distribution the mean does not equal the median.

Try

> median(giveSample(makeExpDist, 0.008, 1000)$le_num)
[1] 87.15622

The median for an $\text{Exponential}(a)$ is $log(2)/a$ while the mean is $1/a$.

> log(2)/0.008
[1] 86.6434
> 1/0.008
[1] 125

Answer 2

This is not the place for debugging R code, so I will focus on the the theory. It seems you are trying to use the quantile function (inverse CDF) of a sample from $\mathsf{Unif}(0,1)$ to get a random sample from $\mathsf{Exp}(\mathrm{rate}=\lambda=0.008)$

As you say, the quantile function is $X = -\log(1-U)/\lambda.$ Notice that $U^\prime = 1 - U \sim\mathsf{Exp}(\lambda)$ also, so it is a little more efficient to use $X = -\log(U)/\lambda.$ (But either method works.) Here is simple R code, which gives the mean of a sample of a million observations from $\mathsf{Exp}(\lambda=0.008)$ as $124.93 \pm 0.25,$ as anticipated.

set.seed(2021)
u = runif(10^6)
x = -log(u)/.008
mean(x)
[1] 124.9272     ## aprx E(X) = 125
2*sd(x)/1000
[1] 0.250065     ## aprx 95% margin of simulation error

mean(-log(1-u)/0.008)
[1] 125.1309     ## alternate method

In the histograms below images and objects are the same color.

R code for figure:

frb = rainbow(8)
cutp.u=seq(0,1, by=.2)
par(mfrow=c(1,2))
 hist(u, prob=T, br=cutp.u, col=frb[c(1,2,3,5,6)])
  curve(dunif(x), add=T, lwd=2)
 cutp.x = -log(cutp.u[6:1])/.008; cutp.x[6]=max(x)
 hist(x, prob=T, br=cutp.x, col=frb[c(6,5,3,2,1)])
  curve(dexp(x,.008), add=T, lwd=2)
par(mfrow=c(1,1))