4

我有兴趣构建一个 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
4

3 回答 3

5

有两个高斯比率的 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 的任何泰勒展开式。

于 2016-02-09T04:34:37.757 回答
3

牢记@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)。对于 Cauchygratio(0,0,1,1)案例NA（对于由 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 的另一个问题的时间？）

于 2016-02-09T02:02:47.953 回答
2

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

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)
于 2016-02-09T00:10:38.280 回答