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')