Gaussian Ratio Distribution: Derivatives wrt underlying $\mu$'s and $\sigma^2$s

Question

I'm working with two independent normal distributions $X$ and $Y$, with means $\mu_x$ and $\mu_y$ and variances $\sigma^2_x$ and $\sigma^2_y$.

I'm interested in the distribution of their ratio $Z=X/Y$. Neither $X$ nor $Y$ has a mean of zero, so $Z$ is not distributed as a Cauchy.

I need to find the CDF of $Z$, and then take the derivative of the CDF with respect to $\mu_x$, $\mu_y$, $\sigma^2_x$ and $\sigma^2_y$.

Does anyone know a paper where these have already been calculated? Or how to do this myself?

I found the formula for the CDF in a 1969 paper, but taking these derivatives will definitely be a huge pain. Maybe someone has already done it or knows how to do it easily? I mainly need to know the signs of these derivatives.

This paper also contains an analytically simpler approximation if $Y$ is mostly positive. I can't have that restriction. However, maybe the approximation has the same sign as the true derivative even outside the parameter range?

Answer 1

Some related papers:

Wiki: `http://en.wikipedia.org/wiki/Ratio_distribution

1. http://www.jstatsoft.org/v16/i04/

2. http://link.springer.com/article/10.1007/s00362-012-0429-2

3. http://mrvar.fdv.uni-lj.si/pub/mz/mz1.1/cedilnik.pdf

Answer 2

Consider using a symbolic math package like Mathematica ,if you have a license, or Sage if you don't.

If your just doing numerical work, you might also just consider numerical differentiation.

While tedious, it does look straight forward. That is, all of the functions involved have easy to compute derivatives. You might use numerical differentiation to test your result when you are done to be sure you have the right formula.

Answer 3

This is the sort of problem that is very easy numerically, and less error prone as well. Since you say you only need the signs, I assume that accurate numerical approximations are more than sufficient for your needs. Here is some code with an example of the derivative against $\mu_x$:

pratio <- function(z, mu_x=1.0, mu_y=1.0,var_x=0.2, var_y=0.2) {
    sd_x <- sqrt(var_x)
    sd_y <- sqrt(var_y)

    a <- function(z) {
        sqrt(z*z/var_x+1/var_y)
    }

    b <- function(z) {
        mu_x*z/var_x + mu_y/var_y
    }

    c <- mu_x^2/var_x + mu_y^2/var_y

    d <- function(z) {
        exp((b(z)^2 - c*a(z)^2)/(2*a(z)^2))
    }


    t1 <- (b(z)*d(z)/a(z)^3)
    t2 <- 1.0/(sqrt(2*pi)*sd_x*sd_y)
    t3 <- pnorm(b(z)/a(z)) - pnorm(-b(z)/a(z))
    t4 <- 1.0/(a(z)^2*pi*sd_x*sd_y)
    t5 <- exp(-c/2.0)
    return(t1*t2*t3 + t4*t5)
}

# Integrates to 1, so probably no typos.
print(integrate(pratio, lower=-Inf, upper=Inf))

cdf_ratio <- function(x, mu_x=1.0, mu_y=1.0,var_x=0.2, var_y=0.2) {
    integrate(function(x) {pratio(x, mu_x, mu_y, var_x, var_y)}, 
        lower=-Inf, upper=x, abs.tol=.Machine$double.eps)$value
} 

# Numerical differentiation here is very easy:
derv_mu_x <- function(x, mu_x=1.0, mu_y=1.0,var_x=0.2, var_y=0.2) {
    eps <- sqrt(.Machine$double.eps)
    left <- cdf_ratio(x, mu_x+eps, mu_y, var_x, var_y)
    right <- cdf_ratio(x, mu_x-eps, mu_y, var_x, var_y)
    return((left - right)/(2*eps))
}