windows(record=T)
n <- 10000
xi <- 1
lambda <- 2
uppertail <- 7
lowertail <- 0.01
gt <- function(alpha,xi,lambda) {
  (1/gamma(alpha)^lambda)*(xi^(alpha-1))
  }
q <- function(alpha) gt(alpha,xi,lambda)
k <- integrate(q,0,10)$value
g <- function(alpha) q(alpha)/k
ga <- function(alpha) alpha*g(alpha)
m <- integrate(ga,0,10)$value
ga2 <- function(alpha) (alpha-m)^2*g(alpha)
v <- integrate(ga2,0,10)$value
mtry <- 1.5
stry <- 0.75
nd <- function(x) dnorm(x,mtry,stry)
betatry <- stry^2/mtry
alphatry <- mtry/betatry
p <- function(x) dgamma(x,shape=alphatry,scale=betatry)
curve(q,0.01,4,xlab="",ylim=c(0.01,1.3),ylab="")
par(new=T)
curve(nd,0.01,4,xlab="",ylim=c(0.01,1.3),ylab="",lty=2)
par(new=T)
curve(p,0.01,4,xlab="",ylim=c(0.01,1.3),ylab="",lty=3)
legend(2.5,1.1,c("Conjugate gamma","Normal","Gamma"),lty=c(1:3))
cat("Ratio of un-normalized density of lower tail at x =",lowertail,
  "with that of normal \n  of mean",mtry,"and s.d.",stry,"is",
  q(lowertail)/nd(lowertail),"\n")
cat("Ratio of un-normalized density of lower tail at x =",lowertail,
  "with that of normal \n  of mean",mtry,"and s.d.",stry,"is",
  q(lowertail)/p(lowertail),"\n")
cat("Ratio of un-normalized density of upper tail at x =",uppertail,
  "with that of normal \n  of mean",mtry,"and s.d.",stry,"is",
  q(uppertail)/nd(uppertail),"\n")
cat("Ratio of un-normalized density of upper tail at x =",uppertail,
  "with that of gamma \n  of same mean and s.d.",
  q(uppertail)/p(uppertail),"\n")
x <- rgamma(n,shape=alphatry,scale=betatry)
kk <- q(x)/p(x)
mm <- (x/k)*q(x)/p(x)
vv <- ((x-mean(mm))^2/k)*q(x)/p(x)
cat("\nImportance sampling\n")
cat("Constant of proportionality",mean(kk),"; true value",k,"\n")
cat("Mean",mean(mm),"; true value",m,"\n")
cat("Variance",mean(vv),"; true value",v,"\n")
w <- q(x)/p(x)
pi <- w/sum(w)
samp <- sample(x,size=n,prob=pi,replace=T)
qmed <- function(x) integrate(g,0.01,x)$value-0.5
med <- uniroot(qmed,lower=0.01,upper=4)$root
cat("\nSampling importance resampling\n")
cat("Mean",mean(samp),"; true value",m,"\n")
cat("Variance",var(samp),"; true value",v,"\n")
cat("Median",median(samp),"; true value",med,"\n")
curve(g,0.01,4,xlab="",ylim=c(0.01,0.6),ylab="")
par(new=T)
plot(density(samp),xlim=c(0.01,4),xlab="",
  ylim=c(0.01,0.6),ylab="",main="")
# for HDR
alpha <- 0.1
le <- (1-alpha)*n
lo <- 1:(n-le)
hi <- (le+1):n
y <- sort(samp)
r <- y[hi]-y[lo]
rm <- min(r)
lom <- min(lo[r==rm])
him <- min(hi[r==rm])
abline(v=y[lom])
abline(v=y[him])