Rejection sampling is quite simple here (see the accepted answer or the end of this post for more efficient alternatives, though):
- Build a rectangle around the density (which is easy here due to bounded support) - the green box in the figure.
- Draw uniform variates in the box (so just two uniform vectors with the right suppport in both directions) - all the dots in the figure, both blank and red.
- Keep those points which are below the density of the double triangular (the red dots in the figure).
- Take the "$x$-axis" coordinates of these points as realizations of the double triangular.
In a picture:
Created with
a <- 0
b <- 2
c <- 1
lambda <- .5
x <- seq(a,b,by=.01)
ddoubletriangular <- function(x,a,b,c,lambda){
p <- lambda/2*(c-a)^2
mu <- 2*(1-p)/(b-c)^2
return(lambda*(x-a)*(x<=c)+mu*(b-x)*(x>c))
}
plot(x,ddoubletriangular(x,a,b,c,lambda),type="l")
upperbound <- max(ddoubletriangular(x,a,b,c,lambda))
lines(c(b,b,a,a,b),c(0,upperbound,upperbound,0,0),col="green",lwd=4)
px <- runif(100,min=a,max=b)
py <- runif(100,min=0,max=upperbound)
points(px,py)
ikeep <- (py < ddoubletriangular(px,a,b,c,lambda))
px <- px[ikeep]
py <- py[ikeep]
points(px,py,col="red",pch=19)
Let's try more than just $N=100$ draws to study the quality of the algorithm:
N <- 1500000
px <- runif(N,min=a,max=b)
py <- runif(N,min=0,max=upperbound)
ikeep <- (py < ddoubletriangular(px,a,b,c,lambda))
px <- px[ikeep]
truehist(px,main="Rejection method",sub=paste("n=",length(px),"of",N))
lines(x,ddoubletriangular(x,a,b,c,lambda),type="l")
This results in
UPDATE:
As forcefully documented in the comments, the above algorithm is, being "off-the shelf", unnecessarily inefficient. Here is an illustration of what @whuber proposes (or so I hope):
First improvement:
Build two rectangular boxes and perform rejection sampling in each region separately (with probability $p$ for the first segment, and $1-p$ for the second), see the next figure. This by construction improves the acceptance probability to 1/2.
Second improvement:
Recognize that the regions in the boxes each are what @whuber calls "rectangular symmetric": the triangle above a segment of the linear density is a rotated version of the acceptance triangle below the density. Hence, we may use the "rejected" samples of the basic algorithm and use these as legitimate realizations for the other end of the segment.
For example, a rejected draw with $x$-coordinate $c+\epsilon$ may be used as an accepted draw for $b-\epsilon$, because the "height" of the rejection region at $c+\epsilon$ is by construction identical to the "height" of the acceptance region at $b-\epsilon$. This has the nice feature of increasing the acceptance rate to 100%.
Here is a way to code that:
a <- 0
b <- 2
c <- 1
lambda <- .5
x <- seq(a,b,by=.002)
p <- lambda/2*(c-a)^2
mu <- 2*(1-p)/(b-c)^2
ddoubletriangular <- lambda*(x-a)*(x<=c)+mu*(b-x)*(x>c)
upperbound1 <- lambda*c
upperbound2 <- mu*x[min(which(x>c))+1]
N <- 1500000
px1 <- runif(N*p,min=a,max=c)
py1 <- runif(N*p,min=0,max=upperbound1)
px2 <- runif(N*(1-p),min=c,max=b)
py2 <- runif(N*(1-p),min=0,max=upperbound2)
ikeep1 <- (py1 < lambda*(px1-a))
px1a <- px1[ikeep1]
px1b <- a+c-px1[!ikeep1]
ikeep2 <- (py2 < mu*(b-px2))
px2a <- px2[ikeep2]
px2b <- c+b-px2[!ikeep2]
px <- c(px1a,px1b,px2a,px2b)
truehist(px,main="Modified rejection method",sub=paste("n=",length(px),"of",N))
lines(x,lambda*(x-a)*(x<=c)+mu*(b-x)*(x>c),type="l")