高斯比率分布:基于和 s 的导数μμσ2σ2

机器算法验证 分布 正态分布 参考 数理统计 累积分布函数
2022-01-26 08:10:25

我正在使用两个独立的正态分布,平均值为以及方差XYμxμyσx2σy2

我对他们的比率的分布感兴趣。的均值都不不是作为柯西分布的。Z=X/YXYZ

我需要找到的 CDF ,然后取 CDF 关于的导数。Zμxμyσx2σy2

有谁知道已经计算过这些的论文?或者自己如何做到这一点?

我在1969 年的一篇论文中找到了 CDF 的公式,但使用这些导数肯定会很痛苦。也许有人已经做过或知道如何轻松做到?我主要需要知道这些衍生物的迹象。

大部分为正,本文还包含一个分析上更简单的近似值。我不能有这个限制。但是,即使在参数范围之外,近似值也可能与真实导数具有相同的符号?Y

3个回答

如果您有许可证,请考虑使用像 Mathematica 这样的符号数学包,如果没有许可证,请考虑使用 Sage。

如果你只是做数值工作,你也可以只考虑数值微分。

虽然乏味,但它确实向前看。也就是说,所有涉及的函数都易于计算导数。完成后,您可以使用数值微分来测试您的结果,以确保您有正确的公式。

这是一种在数值上非常简单的问题,而且也不太容易出错。由于您说您只需要符号,因此我认为准确的数值近似值足以满足您的需求。这是一些代码,其中包含针对的导数示例: μx

pratio <- function(z, mu_x=1.0, mu_y=1.0,var_x=0.2, var_y=0.2) {
    sd_x <- sqrt(var_x)
    sd_y <- sqrt(var_y)

    a <- function(z) {
        sqrt(z*z/var_x+1/var_y)
    }

    b <- function(z) {
        mu_x*z/var_x + mu_y/var_y
    }

    c <- mu_x^2/var_x + mu_y^2/var_y

    d <- function(z) {
        exp((b(z)^2 - c*a(z)^2)/(2*a(z)^2))
    }


    t1 <- (b(z)*d(z)/a(z)^3)
    t2 <- 1.0/(sqrt(2*pi)*sd_x*sd_y)
    t3 <- pnorm(b(z)/a(z)) - pnorm(-b(z)/a(z))
    t4 <- 1.0/(a(z)^2*pi*sd_x*sd_y)
    t5 <- exp(-c/2.0)
    return(t1*t2*t3 + t4*t5)
}

# Integrates to 1, so probably no typos.
print(integrate(pratio, lower=-Inf, upper=Inf))

cdf_ratio <- function(x, mu_x=1.0, mu_y=1.0,var_x=0.2, var_y=0.2) {
    integrate(function(x) {pratio(x, mu_x, mu_y, var_x, var_y)}, 
        lower=-Inf, upper=x, abs.tol=.Machine$double.eps)$value
} 

# Numerical differentiation here is very easy:
derv_mu_x <- function(x, mu_x=1.0, mu_y=1.0,var_x=0.2, var_y=0.2) {
    eps <- sqrt(.Machine$double.eps)
    left <- cdf_ratio(x, mu_x+eps, mu_y, var_x, var_y)
    right <- cdf_ratio(x, mu_x-eps, mu_y, var_x, var_y)
    return((left - right)/(2*eps))
}