windows(record=T)
N <- 5000
mu <- 1
beta <- 1.3348
pump <- c(0.0626,0.1181,0.0937,0.1176,0.6115,
          0.6130,0.8664,0.8661,1.4958,1.9416)
k <- length(pump)
logf <- function(alpha) {
  (alpha-1)*sum(log(pump))-k*log(gamma(alpha))-
  k*alpha*log(beta)-alpha/mu
  }
dlogf <- function(alpha) {
  sum(log(pump))-k*digamma(alpha)-
  k*log(beta)-1/mu
  }
scale <- 4/3
f <- function(alpha) exp(logf(alpha))
n2 <- uniroot(dlogf,lower=0.01,upper=100)$root
h2 <- logf(n2)
n1 <- n2/scale
h1 <- logf(n1)
b1 <- dlogf(n1)
a1 <- h1-n1*b1
n3 <- n2*scale
h3 <- logf(n3)
b3 <- dlogf(n3)
a3 <- h3-n3*b3
#  First plot; Log density
curve(logf,0.01,2,ylim=c(-8,1),xlab="Log density")
abline(h=h2)
abline(a1,b1)
abline(a3,b3)
text(0,h2+0.2,"h2")
text(0.035,a1,"a1+b1x")
text(1.5,a3+b3*1.5,"a3+b3x")
#  End of first plot
d <- exp(h2)
f1 <- function(x) exp(a1+b1*x)
f3 <- function(x) exp(a3+b3*x)
t1 <- (h2-a1)/b1
t3 <- (h2-a3)/b3
#  Second plot; Un-normalized density
curve(f,0.01,2,ylim=c(0,0.1),xlab="Un-normalized Density")
abline(h=d)
par(new=T)
curve(f1,0.01,2,ylim=c(0,0.1),ylab="",xlab="")
par(new=T)
curve(f3,0.01,2,ylim=c(0,0.1),ylab="",xlab="")
abline(v=t1)
abline(v=t3)
text(t1-0.03,0,"t1")
text(t3-0.03,0,"t3")
text(0,d+0.002,"d")
text(0.25,f1(0.25)+0.003,"f1")
text(1.25,f3(1.25)+0.003,"f3")
#  End of second plot
q1 <- exp(a1)*(exp(b1*t1)-1)/b1
q2 <- exp(h2)*(t3-t1)
q3 <- -exp(a3+b3*t3)/b3
const <- q1 + q2 + q3
p <- c(q1,q2,q3)/const
case <- sample(1:3,N,replace=T,prob=p)
v <- runif(N)
w <- (case==1)*(1/b1)*log(q1*b1*exp(-a1)*v+1)+
     (case==2)*(t1+v*(t3-t1))+
     (case==3)*(1/b3)*(log(q3*b3*(v-1))-a3)
dq <- d/const
f1q <- function(x) f1(x)/const
f3q <- function(x) f3(x)/const
h <- function(x) {
       dq +
       (x<t1)*(f1q(x)-dq) +
       (x>t3)*(f3q(x)-dq)
     }
# Third plot; Empirical density
plot(density(w),xlim=c(0.01,2),ylim=c(0,2.5),
  xlab="Empirical density before rejection",
  main="")
par(new=T)
curve(h,0.01,2,ylim=c(0,2.5),ylab="",xlab="")
#  End of third plot
u <- runif(N)
alpha <- w
alpha[u>f(alpha)/(const*h(alpha))] <- NA
alpha <- alpha[!is.na(alpha)]
cat("Acceptances",100*length(alpha)/length(w),"%\n")
int <- integrate(f,0.001,100)$value
fnorm <- function(x) f(x)/int
#  Fourth plot
plot(density(alpha),xlim=c(0.01,2),ylim=c(0,2.5),
  xlab="Empirical density of alpha",main="")
par(new=T)
curve(fnorm,0.01,2,ylim=c(0,2.5),ylab="",xlab="")
#  End of fourth plot