Estimation de l'écart-type d'un ratio en utilisant l'expansion de Taylor

Question

Je suis intéressé de construire une fonction R que je peux utiliser pour tester les limites de l'approximation de Taylor. Je suis conscient qu'il y a des limites à ce que je fais, mais c'est exactement ces limites que je souhaite étudier. J'ai deux variables aléatoires normalement distribuées x et y. x a une moyenne de 7 et un écart-type (SD) de 1. y a une moyenne de 5 et un écart-type de 4.

me.x <- 4; sd.x <- 1 
me.y <- 5; sd.y <- 4 

je sais pour estimer le rapport moyen de y/x, comme celui-

# E(y/x) = E(y)/E(x) - Cov(y,x)/E(x)^2 + Var(x)*E(y)/E(x)^3 
me.y/me.x - 0/me.x^2 + sd.x*me.y/me.x^3 
[1] 1.328125 

Je suis toutefois bloqué sur la façon d'estimer l'écart type du ratio? Je me rends compte que je dois utiliser une expansion de Taylor, mais pas comment l'utiliser.

Faire une simple simulation je reçois

x <- rnorm(10^4, mean = 4, sd = 1); y <- rnorm(10^4, mean = 5, sd = 4) 
sd(y/x) 
[1] 2.027593 
mean(y/x)[1] 
1.362142 

Answer 1

Ces approximations sont peu susceptibles d'être utiles car la distribution ne peut pas avoir un écart-type fini. Regardez comment il est instable:

set.seed(123) 
n <- 10^6 
X <- rnorm(n, me.x, sd.x) 
Y <- rnorm(n, me.y, sd.y) 

sd(head(Y/X, 10^3)) 
## [1] 1.151261 

sd(head(Y/X, 10^4)) 
## [1] 1.298028 

sd(head(Y/X, 10^5)) 
## [1] 1.527188 

sd(Y/X) 
## [1] 1.863168 

qui contraste avec ce qui se passe quand nous essayons la même chose avec une variable aléatoire normale:

sd(head(Y, 10^3)) 
## [1] 3.928038 

sd(head(Y, 10^4)) 
## [1] 3.986802 

sd(head(Y, 10^5)) 
## [1] 3.984113 

sd(Y) 
## [1] 3.999024 

Note: Si vous étiez dans une situation différente, par exemple le dénominateur a un support compact, alors vous pouvez le faire:

library(car) 

m <- c(x = me.x, y = me.y) 
v <- diag(c(sd.x, sd.y)^2) 
deltaMethod(m, "y/x", v) 

Answer 2

Avec les précautions suggérées par @ G.Grothendieck à l'esprit: un mnémotechnique utile pour des produits et des quotients de indépendants X et Y des variables est

CV^2(X/Y) = CV^2(X*Y) = CV^2(X) + CV^2(Y) 

où CV est le coefficient de variation (sd(X)/mean(X)), de sorte CV^2 est Var/mean^2. En d'autres termes

Var(Y/X)/(m(Y/X))^2 = Var(X)/m(X)^2 + Var(Y)/m(Y)^2 

ou réorganisés

sd(Y/X) = sqrt[ Var(X)*m(Y/X)^2/m(X)^2 + Var(Y)*m(Y/X)^2/m(Y)^2 ] 

Pour les variables aléatoires avec la moyenne loin de zéro, c'est une approximation raisonnable.

set.seed(101) 
y <- rnorm(1000,mean=5) 
x <- rnorm(1000,mean=10) 
myx <- mean(y/x) 
sqrt(var(x)*myx^2/mean(x)^2 + var(y)*myx^2/mean(y)^2) ## 0.110412 
sd(y/x) ## 0.1122373 

En utilisant votre exemple est bien pire parce que le CV de Y est proche de 1 - J'ai d'abord pensé qu'il avait l'air correct, mais maintenant je vois qu'il est partial et ne pas capturer la variabilité très bien (je m de brancher aussi dans les valeurs attendues de la moyenne et SD plutôt que leurs valeurs simulées, mais pour un si grand échantillon qui devrait être une partie mineure de l'erreur.)

me.x <- 4; sd.x <- 1 
me.y <- 5; sd.y <- 4 
myx <- me.y/me.x - 0/me.x^2 + sd.x*me.y/me.x^3 
x <- rnorm(1e4,me.x,sd.x); y <- rnorm(1e4,me.y,sd.y) 
c(myx,mean(y/x)) 
sdyx <- sqrt(sd.x^2*myx^2/me.x^2 + sd.y^2*myx^2/me.y^2) 
c(sdyx,sd(y/x))  
## 1.113172 1.197855 

rvals <- replicate(1000, 
    sd(rnorm(1e4,me.y,sd.y)/rnorm(1e4,me.x,sd.x))) 
hist(log(rvals),col="gray",breaks=100) 
abline(v=log(sdyx),col="red",lwd=2) 
min(rvals) ## 1.182698 

Tous les mis en conserve méthode delta approches pour le calcul de la varia nce de Y/X utilise l'estimation ponctuelle pour Y/X (c.-à-d.m (Y/X) = mY/mX), plutôt que l'approximation de second ordre que vous avez utilisée ci-dessus. La construction de formes d'ordre supérieur pour la moyenne et la variance devrait être simple si peut-être fastidieux (un système d'algèbre informatique pourrait aider ...)

mvec <- c(x = me.x, y = me.y) 
V <- diag(c(sd.x, sd.y)^2) 
car::deltaMethod(mvec, "y/x", V) 
##  Estimate  SE 
## y/x  1.25 1.047691 

library(emdbook) 
sqrt(deltavar(y/x,meanval=mvec,Sigma=V)) ## 1.047691 

sqrt(sd.x^2*(me.y/me.x)^2/me.x^2 + sd.y^2*(me.y/me.x)^2/me.y^2) ## 1.047691 

Pour ce que ça vaut, je pris le code dans la réponse de @ SeverinPappadeux et fait en une fonction gratio(mx,my,sx,sy). Pour le cas de Cauchy (gratio(0,0,1,1)), il devient confus et rapporte une moyenne de 0 (qui devrait être NA/divergent) mais signale correctement la variance/std dev comme divergent. Pour les paramètres spécifiés par l'OP (gratio(5,4,4,1)) il donne la moyenne = 1.352176, sd = NA comme ci-dessus. Pour les premiers paramètres que j'ai essayés ci-dessus (gratio(10,5,1,1)) il donne la moyenne = 0.5051581, sd = 0.1141726.

Ces expériences numériques suggèrent fortement pour moi que le rapport de gaussiennes parfois a une variance bien définie, mais je ne sais pas quand (temps pour une autre question sur les mathématiques StackOverflow ou CrossValidated?)

Answer 3

Il est une expression analytique pour le PDF du rapport de deux gaussiennes, faite par David Hinkley (voir par exemple Wikipedia). Nous avons donc pu calculer tous les moments, moyens, etc. Je l'ai tapé et apparemment il n'a clairement pas de second moment fini, donc il n'a pas d'écart-type fini. Remarque, j'ai noté votre Y gaussien comme mon X, et votre X comme mon Y (les formules supposent X/Y). J'ai une valeur moyenne de ratio assez proche de ce que vous avez de la simulation, mais la dernière intégrale est infinie, désolé. Vous pouvez goûter de plus en plus de valeurs, mais de l'échantillonnage std.dev se développe également, comme l'a noté @ G.Grothendieck

library(ggplot2) 

m.x <- 5; s.x <- 4 
m.y <- 4; s.y <- 1 

a <- function(x) { 
    sqrt((x/s.x)^2 + (1.0/s.y)^2) 
} 

b <- function(x) { 
    (m.x*x)/s.x^2 + m.y/s.y^2 
} 

c <- (m.x/s.x)^2 + (m.y/s.y)^2 

d <- function(x) { 
    u <- b(x)^2 - c*a(x)^2 
    l <- 2.0*a(x)^2 
    exp(u/l) 
} 

# PDF for the ratio of the two different gaussians 
PDF <- function(x) { 
    r <- b(x)/a(x) 
    q <- pnorm(r) - pnorm(-r) 

    (r*d(x)/a(x)^2) * (1.0/(sqrt(2.0*pi)*s.x*s.y)) * q + exp(-0.5*c)/(pi*s.x*s.y*a(x)^2) 
} 

# normalization 
nn <- integrate(PDF, -Inf, Inf) 
nn <- nn[["value"]] 

# plot PDF 
p <- ggplot(data = data.frame(x = 0), mapping = aes(x = x)) 
p <- p + stat_function(fun = function(x) PDF(x)/nn) + xlim(-2.0, 6.0) 
print(p) 

# first momentum 
m1 <- integrate(function(x) x*PDF(x), -Inf, Inf) 
m1 <- m1[["value"]] 

# mean 
print(m1/nn) 

# some sampling 
set.seed(32345) 
n <- 10^7L 
x <- rnorm(n, mean = m.x, sd = s.x); y <- rnorm(n, mean = m.y, sd = s.y) 
print(mean(x/y)) 
print(sd(x/y)) 

# second momentum - Infinite! 
m2 <- integrate(function(x) x*x*PDF(x), -Inf, Inf) 

Ainsi, il est impossible de tester toute expansion Taylor pour std.dev.