Eric Fail February 2016

Estimating the Standard Deviation of a ratio using Taylor expansion

I am interested to build a R function that I can use to test the limits of the Taylor series approximation. I am aware that there is limits to what I am doing, but it's exactly those limits I wish to investigate.

I have two normally distributed random variables x and y. x has a mean of 7 and a standard deviation (sd) of 1. y has a mean of 5 and a sd of 4.

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

I know how to estimate the mean ratio of y/x, like this

# 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

I am however stuck on how to estimate the Standard Deviation of the ratio? I realize I have to use a Taylor expansion, but not how to use it.

Doing a simple simulation I get

 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

Answers


G. Grothendieck February 2016

Such approximations are unlikely to be useful since the distribution may not have a finite standard deviation. Look at how unstable it is:

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

Contrast that with what happens when we try the same thing with a normal random variable:

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: If you were in a different situation, e.g. the denominator has compact support, then you could do this:

library(car)

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


Ben Bolker February 2016

With the cautions suggested by @G.Grothendieck in mind: a useful mnemonic for products and quotients of independent X and Y variables is

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

where CV is the coefficient of variation (sd(X)/mean(X)), so CV^2 is Var/mean^2. In other words

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

or rearranging

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

For random variables with the mean well away from zero, this is a reasonable approximation.

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

Using your example is considerably worse because the CV of Y is close to 1 -- I initially thought it looked OK, but now I see that it's biased as well as not capturing the variability very well (I'm also plugging in the expected values of the mean and SD rather than their simulated values, but for such a large sample that should be a minor part of the error.)

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

enter image description here

All the canned delta-method approaches to computing the variance of Y/X use the point estimate for Y/X (i.e. m(Y


Severin Pappadeux February 2016

There is an analytical expression for the PDF of the ratio of two gaussians, done by David Hinkley (e.g. see Wikipedia). So we could compute all momentums, means etc. I typed it and apparently it clearly doesn't have finite second momentum, thus it doesn't have finite standard deviation. Note, I've denoted your Y gaussian as my X, and your X as my Y (formulas assume X/Y). I've got mean value of ratio pretty close to the what you've got from simulation, but last integral is infinite, sorry. You could sample more and more values, but from sampling std.dev is growing as well, as noted by @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)

Thus, it is impossible to test any Taylor expansion for std.dev.

Post Status

Asked in February 2016
Viewed 1,328 times
Voted 12
Answered 3 times

Search




Leave an answer