mode.rho <- function(a1,b1, a2,b2) ifelse(a2 >-1, b2/b1* (a1-1)/(a2+1), Inf) E.rho <- function(a1,b1, a2,b2) ifelse(a2 > 1, b2/b1* a1 /(a2-1), Inf) var.rho <- function(a1,b1, a2,b2) ifelse(a2 > 2, (b2/b1)^2 * ( a1 /(a2-1) * ((a1+1)/(a2-2) - a1/(a2-1))), Inf) sigma.rho <- function(a1,b1, a2,b2) ifelse(a2 > 2, sqrt(var.rho(a1,b1, a2,b2)), Inf) f.rho <- function(rho, a1,b1, a2,b2) { lf = ( a1*log(b1) + a2*log(b2) + (a1-1)*log(rho) + (-a1-a2)*log(b2+rho*b1) - lbeta(a1,a2) ) return(exp(lf)) } x1 = 1; T1 = 1 x2 = 2; T2 = 2 a1 = x1+1; b1 = T1 a2 = x2+1; b2 = T2 rho.max = 8 n = 10^6 cat(sprintf("x1,T1 = %.2f, %.2f; ", x1, T1 )) cat(sprintf("x2,T2 = %.2f, %.2f; \n", x2, T2 )) cat(sprintf("alpha1,beta1 = %d, %d; ", a1, b1 )) cat(sprintf("alpha2,beta2 = %d, %d \n", a2, b2 )) Erho <- E.rho(a1,b1, a2,b2) Srho <- sigma.rho(a1,b1, a2,b2) cat(sprintf("mode = %.3f; E() = %.3f ; sigma %.3f\n", mode.rho(a1,b1, a2,b2), Erho, Srho )) x1.r <- rgamma(n, a1, b1) x2.r <- rgamma(n, a2, b2) rho.r <- x1.r/x2.r Mrho <- mean(rho.r) SDrho <- sd(rho.r) cat(sprintf("MC: mean = %.3f ; sigma %.3f\n", Mrho, SDrho )) rho.r = rho.r[rho.r<rho.max] # only for histogram!! # Warning!! It changes normalization! norma = length(rho.r)/n h <- hist(rho.r, nc=150, plot=FALSE) h$density <- h$density * norma h$counts <- h$counts * norma plot(h, col='cyan', freq=FALSE, main=”, xlim=c(0,rho.max), ylim=c(0,0.52), xlab=expression(rho), ylab=expression(paste('f(',rho,')'))) rho = seq(0, rho.max, len=101) points(rho, f.rho(rho, a1,b1, a2,b2), ty='l', col='blue') text(6,0.46,bquote(x[1] == .(x1) ~ "," ~ T[1] == .(T1) ), cex=1.5, col='blue') text(6,0.41,bquote(x[2] == .(x2) ~ "," ~ T[2] == .(T2) ), cex=1.5, col='blue') Erho.s <- round(Erho, 2) Srho.s <- round(Srho, 2) mode.s <- round(mode.rho(a1,b1, a2,b2),2) Mrho.s <- round(Mrho,2) SDrho.s <- round(SDrho,2) text(6,0.35,bquote(E(rho) == .(Erho.s) ~ "," ~ sigma(rho) == .(Srho.s) ), cex=1.5, col='blue') text(6,0.28,bquote(mode(rho) == .(mode.s) ), cex=1.5, col='red') text(6,0.21,bquote("mean" == .(Mrho.s) ~ "," ~ "std" == .(SDrho.s) ), cex=1.5, col='blue') abline(v=Erho, col='blue') abline(v=mode.rho(a1,b1, a2,b2), col='red',)