使用泰勒展开估计比率的标准偏差

Question

我有兴趣构建一个R函数,我可以用它来测试泰勒级数近似的极限.我知道我正在做的事情是有限的,但这正是我希望调查的那些限制.

我有两个正态分布的随机变量x和y. x平均值为7,标准偏差(sd)为1. y平均值为5,sd值为4.

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

我知道如何估算平均比例y/x,就像这样

# 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

然而,我仍然坚持如何估计比率的标准差？我意识到我必须使用泰勒扩展,但不是如何使用它.

我做了一个简单的模拟

 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

由David Hinkley完成的两位高斯比率的PDF分析表达式(例如参见维基百科).因此我们可以计算所有动量,均值等等.我输入它并且显然它没有有限的第二动量,因此它没有有限的标准偏差.注意,我把你的Y gaussian表示为我的X,你的X表示为我的Y(公式假定为X/Y).我的比率的平均值非常接近你从模拟得到的,但最后的积分是无限的,抱歉.您可以采样越来越多的值,但是从采样std.dev也在增长,正如@ 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)

因此,不可能测试std.dev的任何泰勒展开.