B.3 - Ratio of rates: exact evaluations vs simulation

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