1
$\begingroup$

See edit below as question has been answered.

I want to sample from an exponential distribution with parameter $\lambda>0$ truncated in the tail between $a>0$ and $b>0$, such that $b>a$, where the PDF is flat, due to numerical precision of R, and therefore forces me to sample from the logarithmic transform.

After a logarithmic transform, the truncated "log-PDF" is defined by $\log\lambda-\lambda x$ with $x\in[\log(a),\log(b)]$.

Using (acceptance-)rejection sampling, I can draw $N$ samples (in R) from this log-pdf by sampling two uniform random variables $u_{1}\sim U(\log(a),\log(b))$ and $u_{2}\sim U(\log\lambda-\lambda a,\log(\lambda-\lambda b))$ and accepting $u_{1}$ if the pair $(u_{1},u_{2})$ is located below the log-pdf.

The log-pdf at $u_{1}$ is defined by $f_{\log}(u_{1})=\log\lambda-\lambda a-\frac{\lambda b-\lambda a}{\log(b)-\log(a)}(u_{1}-\log(a))$, therefore if $u_{2}\leq f_{\log}(u_{1})$, then $u_{1}$ is accepted.

My R implementation is as follows for $N = 1000$, $\lambda=5$, $a=50$ and $b=55$, such that a comparison can be made by directly drawing from an exponential distribution(; however, the goal is to be able to set, e.g., $a=180$ and $b=200$, where the exponential distribution is flat in R):

N <- 10000     # number of samples
M <- 1000      # number of plot points
lambda <- 5    # parameter
a <- 50        # left truncation point
b <- 55        # right truncation point
x <- rep(0,N)  # empty vector
i <- 1         # count variable
while(i <= N){ # loop until count hits number of samples
  u1 <- runif(1,log(a),log(b))                                              # uniform random variable 1
  u2 <- runif(1,log(lambda)-lambda*b,log(lambda)-lambda*a)                  # uniform random variable 2
  f <- log(lambda)-lambda*a-(lambda*b-lambda*a)/(log(b)-log(a))*(u1-log(a)) # log pdf at u1
  if(u2 <= f){ # check if u1 is accepted
    x[i] <- u1 # save u1
    i <- i+1   # update count variable
  }
}
hist(exp(x))                                                    # plot histogram after applying exponent
plot(seq(a,b,length=M),dexp(seq(a,b,length=M),lambda),type='l') # plot exponential distribution

Unfortunately, the results significantly deviate and therefore my question(s) is/are whether my understanding of using a logarithmic transform and/or rejection sampling is (in)correct (or whether there is an error in my R code)? enter image description hereenter image description here

Edit

The second uniform random variable $u_{2}$ should also undergo a logarithmic transformation, such that the acceptance ratio is correct and thus $u_{1}$ is accepted as follows (reposted the while loop)

 while(i <= N){
  u1 <- runif(1,log(a),log(b))
  u2 <- log(runif(1L))
  fu1 <- log(lambda)-lambda*a-(lambda*b-lambda*a)/(log(b)-log(a))*(u1-log(a))
  floga <- log(lambda)-lambda*a-(lambda*b-lambda*a)/(log(b)-log(a))*(log(a)-log(a)) # log pdf at log(a)
  if(u2 <= fu1-floga){
    x[i] <- u1
    i <- i+1
  }
}
$\endgroup$

0

Browse other questions tagged or ask your own question.