# Simulating ratio of two independent normal variables

I have two non centered independent normally distributed variables $X$ and $Y$. Specifically $$X \sim N(\mu_1, \sigma_1)$$ $$Y \sim N(\mu_2, \sigma_2)$$

I am interested in sampling from the ratio $X/Y$. I know that the ratio is Cauchy distributed if the denominator distribution is centered at 0 which mine is not. My question is is there an approximate distribution that I can sample the ratio $X/Y$ from using the means and standard deviations from the normal distributions?

– whuber
Jul 11, 2018 at 22:26
• If you just want to simulate the ratio, why not simply do that? i.e. simulate X and Y and compute X/Y? Trying to compute the distribution and simulate that seems like a lot of work by comparison Jul 11, 2018 at 23:21

If you are only interested in sampling from this ratio distribution, you can do that directly without knowing the distribution function of the ratio (exactly or approximately). Here is how I would do it in R:

rmyratio <- function(n, mean, sd) {
rnorm(n, mean[1], sd[1]) / rnorm(n, mean[2], sd[2])
}
X  <-  rmyratio(1000, 0:1, c(1,5))
hist(X)


If you for some other purpose want to know the exact distribution, it can be found in https://www.amazon.com/Probability-Distributions-Involving-Gaussian-Variables/dp/0387346570/ref=sr_1_1?s=books&ie=UTF8&qid=1531341533&sr=1-1&keywords=gaussian+distribution+engineering (I am sure, but didn't look)

There is plenty of information about the ratio you are looking for at this nice Wikipedia link.

However, if you just want to sample, you can do that as well, easily using R code...

# set up variables
num.samples <- 100000
mu1 <- 1
mu2 <- 2
sd1 <- 4
sd2 <- 9

# simulate norms
norms1 <- rnorm(n=num.samples, mean=mu1, sd=sd1)
norms2 <- rnorm(n=num.samples, mean=mu2, sd=sd2)

# do ratio
ratio.dist <- norms1 / norms2

# to sample, use sample()
sample(ratio.dist, size=1)

# You can see the density plot as well...
plot(density(ratio.dist))

• I added my contribution after @kjetil, so you should accept their answer if it is what you are looking for.
– ERT
Jul 11, 2018 at 20:45

An exact expression for the ratio of two correlated normal random variables was given by Hinkley in 1969. The notation is taken from the stackexchange post

For $$Z = \frac{X}{Y}$$ with $$\begin{bmatrix}X\\Y\end{bmatrix} \sim N\left(\begin{bmatrix} \mu_x \\ \mu_y \end{bmatrix} , \begin{bmatrix} \sigma_x^2 & \rho \sigma_x \sigma_y \\ \rho \sigma_x \sigma_y & \sigma_y^2 \end{bmatrix} \right)$$

The exact probability density function (pdf) result is:

$$f(z) = \frac{b(z)d(z)}{a(z)^3} \frac{1}{\sqrt{2\pi} \sigma_X\sigma_Y} \left[ \Phi \left( \frac{b(z)}{\sqrt{1-\rho^2}a(z)} \right) - \Phi \left( - \frac{b(z)}{\sqrt{1-\rho^2}a(z)} \right) \right] + \frac{\sqrt{1-\rho^2}}{\pi \sigma_X \sigma_Y a(z)^2} \exp \left( -\frac{c}{2(1-\rho^2)}\right)$$ with $$\begin{array}{} a(z) &=& \left( \frac{z^2}{\sigma_X^2} - \frac{2 \rho z}{\sigma_X \sigma_Y} + \frac{1}{\sigma_Y^2} \right) ^{\frac{1}{2}} \\ b(z) &=& \frac{\mu_X z}{ \sigma_X^2} - \frac{\rho (\mu_X+ \mu_Y z)}{ \sigma_X \sigma_Y} + \frac{\mu_Y}{\sigma_Y^2} \\ c &=& \frac{\mu_X^2}{\sigma_Y^2} - \frac{2 \rho \mu_X \mu_Y + }{\sigma_X \sigma_Y} + \frac{\mu_Y^2}{\sigma_Y^2}\\ d(z) &=& \text{exp} \left( \frac {b(z)^2 - c a(z) ^2}{2(1-\rho^2)a(z)^2}\right) \end{array}$$

And an approximation of the CDF based on an asymptotic behaviour is: (for $$\theta_Y/\sigma_Y \to \infty$$): $$F(z) \to \Phi\left( \frac{z - \mu_X/\mu_Y}{\sigma_X \sigma_Y a(z)/\mu_Y} \right)$$

You end up with the Delta method result when you insert the approximation $$a(z) = a(\mu_X/\mu_Y)$$ $$a(z) \sigma_X \sigma_Y /\mu_Y \approx a(\mu_X/\mu_Y) \sigma_X \sigma_Y /\mu_Y = \left( \frac{\mu_X^2\sigma_Y^2}{\mu_Y^4} - \frac{2 \mu_X \sigma_X \sigma_Y}{\mu_Y^3} + \frac{\sigma_X^2}{\mu_Y^2} \right) ^{\frac{1}{2}}$$