1
$\begingroup$

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!

$\endgroup$
4
  • $\begingroup$ Where did the "$ln(.5)$" come from?? The mean of your distribution is $1/a.$ $\endgroup$
    – whuber
    Jul 8, 2021 at 19:37
  • 1
    $\begingroup$ @whuber I think OP just confused 50th percentile (plugging in $p=0.5$ into their quantile function) for the mean. $\endgroup$
    – bdeonovic
    Jul 8, 2021 at 19:49
  • $\begingroup$ @bdeonovic Good interpretation! (+1) $\endgroup$
    – whuber
    Jul 8, 2021 at 20:34
  • $\begingroup$ You're right, I did confuse the two terms. The mean is $1/.008=125$ which is what I was looking for. Used the wrong integration technique for mean. $\endgroup$ Jul 9, 2021 at 20:11

2 Answers 2

1
$\begingroup$

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
$\endgroup$
1
$\begingroup$

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.

enter image description here

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))
$\endgroup$
1
  • $\begingroup$ Thank you for the graphical demonstration. $\endgroup$ Jul 9, 2021 at 20:13

Your Answer

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

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