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

28

Pracuję z dwoma niezależnymi rozkładami normalnymi i Y , ze średnimi μ x i μ y oraz wariancjami σ 2 x i σ 2 y .XYμxμyσx2σy2

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

Muszę znaleźć CDF dla , a następnie wziąć pochodną CDF w odniesieniu do μ x , μ y , σ 2 x i σ 2 y .Zμxμyσx2σy2

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

ABC
źródło
4
Dodałem dla ciebie. Napisałeś „sigma”, ale wspomniałeś, że były to wariancje, więc uczyniłem je kwadratami sigma. Upewnij się, że nadal mówi, o co chcesz zapytać. TEX
gung - Przywróć Monikę
3
en.wikipedia.org/wiki/Ratio_distribution ma funkcję gęstości prawdopodobieństwa.
Douglas Zare
2
To ten sam plik PDF, co w powyższym artykule. Próbuję wziąć pochodną CDF w odniesieniu do bazowej mus i sigma.
ABC
2
Formuła pdf znaleziona przez Davida Hinkleya jest całkowicie zamknięta. Możesz więc brać te instrumenty pochodne, krok po kroku. Tak naprawdę jestem ciekawy sensu robienia takich pochodnych, ponieważ nie ma powodu, dla którego znak powinien być stały równomiernie w stosunku do liczb rzeczywistych ...
Xi'an
2
X/Y

Odpowiedzi:

1

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

Kwant
źródło
5
Witamy na stronie @Quantum. Czy mógłbyś podać krótkie streszczenie tych artykułów, aby czytelnicy mogli ocenić, czy są tym, czego szukają, bez konieczności otwierania i czytania każdego z nich?
gung - Przywróć Monikę
Z=X/Y
4
Kwantowe - to nie rozwiązuje problemu @ gung. Odpowiedzi zawierające tylko łącze zazwyczaj nie są akceptowane. Gung zapytał, czy możesz „podać krótkie streszczenie tych dokumentów” (co oznacza „w swojej odpowiedzi”). Twój zbiorczy opis w komentarzu nie jest wystarczający. Podaj krótki opis każdego linku (jeśli to możliwe, indywidualnie, a nie zbiorowo), który wskaże, dlaczego go umieściłeś / dlaczego jest on istotny. W tej chwili potencjalnie użyteczna odpowiedź może zostać przekształcona w komentarz - tak jak to już miało miejsce w przypadku wcześniejszych odpowiedzi tylko na link na to pytanie.
Glen_b
XYZ=XYxyp(x,y)dxdy
y
0

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łę.

Dave31415
źródło
0

μ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))
} 
AaronDefazio
źródło