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

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

这样的近似值不太可能有用，因为分布可能没有有限的标准偏差。看看它有多不稳定：

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

对比一下，当我们用一个正常的随机变量尝试同样的事情时会发生什么：

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

注意：如果您处于不同的情况，例如分母有紧凑的支持，那么你可以这样做：

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

牢记@G.Grothendieck 建议的注意事项：独立 X 和 Y 变量的乘积和商的有用助记符是

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

其中CV 是变异系数（sd(X)/mean(X)），所以CV^2 是Var/mean^2。换句话说

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

或重新排列

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

对于均值远离零的随机变量，这是一个合理的近似值。

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

使用您的示例要差得多，因为 Y 的 CV 接近 1 - 我最初认为它看起来不错，但现在我发现它有偏差并且没有很好地捕捉可变性（我也在插入平均值和 SD 的期望值而不是它们的模拟值，但对于如此大的样本，应该是误差的一小部分。）

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

所有计算 Y/X 方差的固定 delta 方法都使用 Y/X 的点估计（即 m(Y/X) = mY/mX），而不是您上面使用的二阶近似值.为均值和方差构建高阶形式应该很简单，如果可能很乏味（计算机代数系统可能会有所帮助......）

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

为了它的价值，我把@SeverinPappadeux 的答案中的代码做成了一个函数gratio(mx,my,sx,sy)。对于 Cauchy 案例 (gratio(0,0,1,1))，它会感到困惑并报告平均值为 0（应该是 NA/divergent），但正确地将方差/标准差报告为发散的。对于由 OP (gratio(5,4,4,1)) 指定的参数，它给出 mean=1.352176, sd=NA 如上所述。对于我在上面尝试的第一个参数 (gratio(10,5,1,1))，它给出了 mean=0.5051581, sd=0.1141726。

这些数值实验强烈地向我表明，高斯的比率有时具有明确定义的方差，但我不知道什么时候（关于 Math StackOverflow 或 CrossValidated 的另一个问题的时间？）

Answer 3

有两个高斯比率的PDF的解析表达式，完成 David Hinkley（例如，参见 Wikipedia）。所以我们可以计算所有的动量、平均值等。我输入了它，显然它没有有限的第二动量，因此它没有有限的标准偏差。请注意，我将您的 Y 高斯表示为我的 X，将您的 X 表示为我的 Y（公式假设 X/Y）。我得到的比率平均值非常接近你从模拟中得到的值，但是最后一个积分是无限的，抱歉。您可以采样越来越多的值，但正如@G.Grothendieck 所指出的那样，std.dev 的采样也在增长

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 的任何泰勒展开式。