<s>[INST] 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?
 [/INST] 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}}$$
</s>