n=10^7
x1 = 1
x2 = 1
lambda1 = rgamma(n, x1+1, 1)
lambda2 = rgamma(n, x2+1, 1)
rho = lambda1/lambda2
E.rho = mean(rho)
sigma.rho = sd(rho)

max.rho.hist = 8
rho = rho[rho<max.rho.hist] 
hist(rho, nc=150, col='cyan', freq=FALSE, xlim=c(0,max.rho.hist), main='',
     xlab=expression(paste(rho, ' = ', lambda[1], '/', lambda[2])) )

# dummy histogram for rough evaluation of the mode
h.rx <- hist(rho, nc=1000, plot=FALSE )
mode = h.rx$mids[which.max(h.rx$density)]

abline(v=mode, col='red')
abline(v=E.rho, col='blue')

p.overflow = (n - length(rho))/n * 100
cat(sprintf("fraction of overflows %.2f%%\n", p.overflow))
text(6,0.63-0.05,expression(paste(x[1], ' =  ', x[2], ' = 1')),
     cex=1.6,col='blue')
text(6,0.54-0.05,sprintf("mean = %.2f;  std = %.2f", E.rho, sigma.rho),
     cex=1.5, col='blue')
text(6,0.46-0.05, sprintf("[ Overflow: %.2f%% ]", p.overflow),
     cex=1.5, col='gray')
text(6,0.37-0.05, sprintf("[ Mode = %.2f ]", mode), cex=1.5, col='red')