Rozkład współczynnika Gaussa: Pochodne wrt leżące u podstaw

Pracuję z dwoma niezależnymi rozkładami normalnymi i Y , ze średnimi μ x i μ y oraz wariancjami σ 2 x i σ 2 y .$X$$Y$$\mu_x$$\mu_y$$\sigma^2_x$$\sigma^2_y$

Jestem zainteresowany w dystrybucji ich stosunek . Ani X, ani Y nie ma średniej zero, więc Z nie jest dystrybuowane jako Cauchy.$Z=X/Y$$X$$Y$$Z$

Muszę znaleźć CDF dla , a następnie wziąć pochodną CDF w odniesieniu do μ x , μ y , σ 2 x i σ 2 y .$Z$$\mu_x$$\mu_y$$\sigma^2_x$$\sigma^2_y$

Czy ktoś zna papier, w którym zostały one już obliczone? Lub jak to zrobić sam?

Wzór na CDF znalazłem w artykule z 1969 roku , ale przyjmowanie tych pochodnych z pewnością będzie ogromnym bólem. Może ktoś już to zrobił lub wie, jak to zrobić z łatwością? Potrzebuję głównie znać znaki tych pochodnych.

Ten artykuł zawiera także analitycznie prostsze przybliżenie, jeśli jest w większości dodatnie. Nie mogę mieć tego ograniczenia. Może jednak aproksymacja ma ten sam znak, co prawdziwa pochodna, nawet poza zakresem parametrów?$Y$

Niektóre powiązane dokumenty:

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

Rozważ użycie symbolicznego pakietu matematycznego, takiego jak Mathematica, jeśli masz licencję, lub Sage, jeśli nie masz.

Jeśli po prostu wykonujesz pracę numeryczną, możesz również rozważyć różnicowanie numeryczne.

Choć nużący, wygląda prosto. Oznacza to, że wszystkie zaangażowane funkcje mają łatwe do obliczenia pochodne. Możesz użyć numerycznego różnicowania do przetestowania wyniku, gdy skończysz, aby upewnić się, że masz właściwą formułę.

$\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))
}