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?
 A: 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)
A: 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))

A: 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}}$$
