The distribution of the ratio of two independent normals and a normal approximation

Let \(X\sim N(\mu_x,\sigma_x^2)\) and \(Y\sim N(\mu_y,\sigma_y^2)\) be two independent random variables, and \(Z=X/Y\). Define \(\beta = \mu_x/\mu_y\), \(\rho=\sigma_y/\sigma_x\), and \(\delta = \sigma_y/\mu_y\). Then, the probability density function of the ratio \(Z\) is given by: \[f_Z(z;\beta,\rho,\delta_y) = \frac{\rho}{\pi(1+\rho^2 z^2)}\Biggl\{\exp\left[-\frac{\rho^2\beta^2+1}{2\delta_y^2}\right] + \sqrt{\frac{\pi}{2}} q \text{erf}\left(\frac{q}{\sqrt{2}}\right)\exp\left[-\frac{\rho^2(z-\beta)^2}{2\delta_y^2(1+\rho^2 z^2)}\right] \Biggr\},\] where \[q=\frac{1+\beta\rho^2 z}{\delta_y\sqrt{1+\rho^2z^2}}.\] This density is heavy tailed and has no moments. The shape of this density can be unimodal, bimodal, symmetric, asymmetric, and/or even similar to a normal distribution close to its mode. [1] found that the density \(f_Z\) can be approximated with a normal distribution on an interval of a certain probability, provided that the coeffiencient of variation of \(Y\), \(\delta_y\), is smaller than \(0.1\) (See Theorem 1 of [1] for more precise conditions). [1] proposed a normal approximation to the distribution of the ratio of two independent normals as \(Z\stackrel{\cdot}{\sim}N(\mu_z,\sigma_z^2)\), where \(\mu_z=\beta = \mu_x/\mu_y\) and \(\sigma_z^2 = \delta_y^2 (\rho^{-2}+\beta^2)\). The following R code shows the shapes of the distribution of \(Z\) as well as the accuracy of the normal approximation proposed by [1]. Note that only the first three examples fall into the range \(\delta_y\leq 0.1\).

References.

  1. On the existence of a normal approximation to the distribution of the ratio of two independent normal random variables
rm(list=ls())
# Error function
erf <- Vectorize(function(x) 2*pnorm(sqrt(2)*x) - 1)

# probability density function of Z
fZ <- function(z,beta,rho,deltay){
  q <- (1+beta*rho^2*z)/(deltay*sqrt(1+rho^2*z^2))
  den <- rho/(pi*(1+rho^2*z^2))*( exp(-(rho^2*beta^2+1)/(2*deltay^2))  + 
                                    sqrt(pi/2)*q*erf(q/sqrt(2))*exp(-0.5*(rho^2*(z-beta)^2)/(deltay^2*(1+rho^2*z^2))) )
  return(den)
}

# normal approximation
dnormapp <- function(z,beta,rho,deltay) dnorm(z,beta,deltay*sqrt(rho^(-2)+beta^2))

# Checking that fZ integrates 1 for some parameter values
tempf <- Vectorize(function(z)  fZ(z,1.7,1.3,0.37))
integrate(tempf,-Inf,Inf)
## 1 with absolute error < 8.5e-08
#####################################################################
# Shape of the density for some parameter values
#####################################################################

# Very very similar to a normal distribution
beta0 <- 2; rho0 <- 0.04; deltay0 <- 0.01;
tempf <- Vectorize(function(z)  fZ(z,beta0,rho0,deltay0))
tempa <- Vectorize(function(z) dnormapp(z,beta0,rho0,deltay0))
curve(tempf,1,3,n=1000,xlab="z",ylab="Density",lwd=2)
curve(tempa,1,3,n=1000,xlab="z",ylab="Density",lty=2,add=T,col="blue",lwd=2)
legend(2.5,1.25, c("fZ","Normal"),
       text.col = "black", lty = c(1, 2),
       merge = TRUE, bg = "gray90")

# Very similar to a normal distribution
beta0 <- 2; rho0 <- 0.05; deltay0 <- 0.05;
tempf <- Vectorize(function(z)  fZ(z,beta0,rho0,deltay0))
tempa <- Vectorize(function(z) dnormapp(z,beta0,rho0,deltay0))
curve(tempf,-2,6,n=1000,xlab="z",ylab="Density",lwd=2)
curve(tempa,-2,6,n=1000,xlab="z",ylab="Density",lty=2,add=T,col="blue",lwd=2)
legend(4, 0.3, c("fZ","Normal"),
       text.col = "black", lty = c(1, 2),
       merge = TRUE, bg = "gray90")

# Similar to a normal distribution
beta0 <- 2; rho0 <- 0.05; deltay0 <- 0.1;
tempf <- Vectorize(function(z)  fZ(z,beta0,rho0,deltay0))
tempa <- Vectorize(function(z) dnormapp(z,beta0,rho0,deltay0))
curve(tempf,-5,10,n=1000,xlab="z",ylab="Density",lwd=2)
curve(tempa,-5,10,n=1000,xlab="z",ylab="Density",lty=2,add=T,col="blue",lwd=2)
legend(6, 0.2, c("fZ","Normal"),
       text.col = "black", lty = c(1, 2),
       merge = TRUE, bg = "gray90")

# Starting to depart from normality
beta0 <- 2; rho0 <- 0.1; deltay0 <- 0.25;
tempf <- Vectorize(function(z)  fZ(z,beta0,rho0,deltay0))
tempa <- Vectorize(function(z) dnormapp(z,beta0,rho0,deltay0))
curve(tempf,-7,12,n=1000,xlab="z",ylab="Density",lwd=2)
curve(tempa,-7,12,n=1000,xlab="z",ylab="Density",lty=2,add=T,col="blue",lwd=2)
legend(8, 0.1, c("fZ","Normal"),
       text.col = "black", lty = c(1, 2),
       merge = TRUE, bg = "gray90")

# Clear departure from normality
beta0 <- 2; rho0 <- 0.5; deltay0 <- 0.5;
tempf <- Vectorize(function(z)  fZ(z,beta0,rho0,deltay0))
tempa <- Vectorize(function(z) dnormapp(z,beta0,rho0,deltay0))
curve(tempf,-7,12,n=1000,xlab="z",ylab="Density",lwd=2)
curve(tempa,-7,12,n=1000,xlab="z",ylab="Density",lty=2,add=T,col="blue",lwd=2)
legend(8, 0.2, c("fZ","Normal"),
       text.col = "black", lty = c(1, 2),
       merge = TRUE, bg = "gray90")

# Clear departure from normality
beta0 <- 0.1; rho0 <- 1; deltay0 <- 5;
tempf <- Vectorize(function(z)  fZ(z,beta0,rho0,deltay0))
tempa <- Vectorize(function(z) dnormapp(z,beta0,rho0,deltay0))
curve(tempf,-7,12,n=1000,xlab="z",ylab="Density",lwd=2)
curve(tempa,-7,12,n=1000,xlab="z",ylab="Density",lty=2,add=T,col="blue",lwd=2)
legend(8, 0.2, c("fZ","Normal"),
       text.col = "black", lty = c(1, 2),
       merge = TRUE, bg = "gray90")

# Bimodal!
beta0 <- 12; rho0 <- 2; deltay0 <- 10;
tempf <- Vectorize(function(z)  fZ(z,beta0,rho0,deltay0))
tempa <- Vectorize(function(z) dnormapp(z,beta0,rho0,deltay0))
curve(tempf,-8,12,n=1000,xlab="z",ylab="Density",lwd=2)
curve(tempa,-8,12,n=1000,xlab="z",ylab="Density",lty=2,add=T,col="blue",lwd=2)
legend(8, 0.2, c("fZ","Normal"),
       text.col = "black", lty = c(1, 2),
       merge = TRUE, bg = "gray90")