ratio2normals <- function(x, mean1,mean2,sd1,sd2,rho){
# A function to compute ratio of 2 normals
# R code written by Ravi Varadhan 
# May 25, 2007
# Based on the paper by Pham Gia et al., Comm in Stats (2006)
A <- 1 / (2*pi*sd1*sd2*sqrt(1-rho^2))
exponent.num <- -sd2^2*mean1^2 - sd1^2*mean2^2 + 2*rho*sd1*sd2*mean1*mean2
exponent.denom <- 2*(1-rho^2)*sd1^2*sd2^2
K <- A * exp(exponent.num/exponent.denom)
t2x.num <- -sd2^2*mean1*x - sd1^2*mean2 + rho*sd1*sd2*(mean2*x + mean1)
t2x.denom <- sd1*sd2*sqrt(2*(1-rho^2)*(sd2^2*x^2 - 2*rho*x*sd1*sd2 + sd1^2))
t2x <- t2x.num / t2x.denom
erf.term <- 2 * pnorm(sqrt(2) * t2x) - 1
Ft2x <- sqrt(pi) * t2x * exp(t2x^2) * erf.term + 1
fx <- K * Ft2x * 2 * (1 - rho^2) * sd1^2 * sd2^2 / (sd2^2 * x^2 + sd1^2 - 2*x*rho*sd1*sd2)
return(fx)
}


mean1 <- 75.25
mean2 <- 71.58 
sd1 <- 6.25
sd2 <- 5.45
rho <- 0.76

x <- seq(0.5,1.5, length=100)
y <- ratio2normals(x,mean1,mean2, sd1,sd2,rho)
plot(x,y, type="l")


# compute the mean and variance via quadrature
m1 <- function(x, mean1, mean2, sd1, sd2, rho){
x * ratio2normals(x, mean1, mean2, sd1, sd2, rho)
}

m2 <- function(x, mean1, mean2, sd1, sd2, rho){
x^2 * ratio2normals(x, mean1, mean2, sd1, sd2, rho)
}

m.1 <- integrate(m1, lower=-Inf, upper=Inf, mean1=mean1, mean2=mean2, sd1=sd1, sd2=sd2, rho=rho, rel.tol=1.e-07)$val

m.2 <- integrate(m2, lower=-Inf, upper=Inf, mean1=mean1, mean2=mean2, sd1=sd1, sd2=sd2, rho=rho, rel.tol=1.e-07)$val

mean <- m.1
variance <- m.2 - m.1^2