n = 10^7 rM = 100 r1 = runif(n, 0, rM) r2 = runif(n, 0, rM) rho = r1/r2 rho.h <- rho[rho<5] # for the histogram norma = length(rho.h)/n # normalization h <- hist(rho.h, nc=100, plot=FALSE) h$density <- h$density * norma h$counts <- h$counts * norma plot(h, col='cyan', freq=FALSE, main='', xlim=c(0,5), ylim=c(0,0.52), xlab=expression(rho), ylab=expression(paste('f(',rho,')'))) abline(v=1, col='red') # r1.check = rho*r2 # hist(r1.check, nc=200, col='blue', freq=FALSE, xlim=c(0,rM))