如何找到给定分布的区间概率？答案

【问题标题】：How to find interval prbability for a given distribution?如何找到给定分布的区间概率？
【发布时间】：2016-12-07 21:42:59
【问题描述】：

假设我有一些数据并将它们拟合到gamma 分布，如何找到Pr(1 < x <= 1.5) 的区间概率，其中 x 是样本外数据点？

require(fitdistrplus)

a <- c(2.44121289,1.70292449,0.30550832,0.04332383,1.0553436,0.26912546,0.43590885,0.84514809,
0.36762336,0.94935435,1.30887437,1.08761895,0.66581035,0.83108270,1.7567334,1.00241339,
0.96263021,1.67488277,0.87400413,0.34639636,1.16804671,1.4182144,1.7378907,1.7462686,
1.7427784,0.8377457,0.1428738,0.71473956,0.8458882,0.2140742,0.9663167,0.7933085,
0.0475603,1.8657773,0.18307362,1.13519144)

fit <- fitdist(a, "gamma",lower = c(0, 0))

【问题讨论】：

x 是一个新的、独立的数据点吗？还是 x 是一个参数估计（即您正在寻找后验概率，例如，形状参数在 1 和 1.5 之间）？
@JacobSocolar 不，x 是数据点。

标签： r probability distribution fitdistrplus

【解决方案1】：

有人不喜欢我上面的方法，这是以 MLE 为条件的；现在让我们看一些无条件的东西。如果我们采用直接集成，我们需要三重集成：一个用于shape，一个用于rate，最后一个用于x。这并不吸引人。我将只生成蒙特卡洛估计值。

根据中心极限定理，MLE 呈正态分布。 fitdistrplus::fitdist 不会给出标准错误，但我们可以使用 MASS::fitdistr 来执行精确推断。

fit <- fitdistr(a, "gamma", lower = c(0,0))

b <- fit$estimate
#   shape     rate 
#1.739737 1.816134 

V <- fit$vcov  ## covariance
          shape      rate
shape 0.1423679 0.1486193
rate  0.1486193 0.2078086

现在我们要从参数分布中采样，得到目标概率的样本。

set.seed(0)
## sample from bivariate normal with mean `b` and covariance `V`
## Cholesky method is used here
X <- matrix(rnorm(1000 * 2), 1000)  ## 1000 `N(0, 1)` normal samples
R <- chol(V)  ## upper triangular Cholesky factor of `V`
X <- X %*% R  ## transform X under desired covariance
X <- X + b  ## shift to desired mean
## you can use `cov(X)` to check it is very close to `V`

## now samples for `Pr(1 < x < 1.5)`
p <- pgamma(1.5, X[,1], X[,2]) - pgamma(1, X[,1], X[,2])

我们可以制作p 的直方图（如果需要，还可以进行密度估计）：

hist(p, prob = TRUE)

现在，我们经常需要预测变量的样本均值：

mean(p)
# [1] 0.1906975

【讨论】：

不错！在这种情况下，参数后验分布的多元正态逼近有多安全？ (+1)
不错的方法，但我注意到fitdistrplus::fitdist 通过最大似然参数提供标准误差。

【解决方案2】：

这里有一个示例，它使用 MCMC 技术和贝叶斯推理模式来估计新观测值落在区间 (1:1.5) 中的后验概率。这是一个无条件估计，与通过将伽马分布与最大似然参数估计相结合获得的条件估计相反。

此代码要求在您的计算机上安装 JAGS（免费且易于安装）。

library(rjags)

a <- c(2.44121289,1.70292449,0.30550832,0.04332383,1.0553436,0.26912546,0.43590885,0.84514809,
       0.36762336,0.94935435,1.30887437,1.08761895,0.66581035,0.83108270,1.7567334,1.00241339,
       0.96263021,1.67488277,0.87400413,0.34639636,1.16804671,1.4182144,1.7378907,1.7462686,
       1.7427784,0.8377457,0.1428738,0.71473956,0.8458882,0.2140742,0.9663167,0.7933085,
       0.0475603,1.8657773,0.18307362,1.13519144)

# Specify the model in JAGS language using diffuse priors for shape and scale
sink("GammaModel.txt")
cat("model{

    # Priors
    shape ~ dgamma(.001,.001)
    rate ~ dgamma(.001,.001)

    # Model structure
    for(i in 1:n){
    a[i] ~ dgamma(shape, rate)
    }
    }
    ", fill=TRUE)
sink()

jags.data <- list(a=a, n=length(a))

# Give overdispersed initial values (not important for this simple model, but very important if running complicated models where you need to check convergence by monitoring multiple chains)
inits <- function(){list(shape=runif(1,0,10), rate=runif(1,0,10))}

# Specify which parameters to monitor
params <- c("shape", "rate")

# Set-up for MCMC run
nc <- 1   # number of chains
n.adapt <-1000   # number of adaptation steps
n.burn <- 1000    # number of burn-in steps
n.iter <- 500000  # number of posterior samples
thin <- 10   # thinning of posterior samples

# Running the model
gamma_mod <- jags.model('GammaModel.txt', data = jags.data, inits=inits, n.chains=nc, n.adapt=n.adapt)
update(gamma_mod, n.burn)
gamma_samples <- coda.samples(gamma_mod,params,n.iter=n.iter, thin=thin)
# Summarize the result
summary(gamma_samples)

# Compute improper (non-normalized) probability distribution for x
x <- rep(NA, 50000)
for(i in 1:50000){
  x[i] <- rgamma(1, gamma_samples[[1]][i,1], rate = gamma_samples[[1]][i,2])
}

# Find which values of x fall in the desired range and normalize.
length(which(x>1 & x < 1.5))/length(x)

回答： Pr(1 < x <= 1.5) = 0.194 非常接近条件估计，但不能保证通常是这种情况。

【讨论】：

@ZheyuanLi 没错，矢量化的方式更好。我循环的原因是因为我不确定 R 的 RNG 是否以这种方式矢量化，并且因为我不关心效率（所以我没有检查）。至于 rgamma 与 pgamma，rgamma 方法是一种万能的方法，可以从 MCMC 输出中获取派生参数的完整后验分布（对于任何任意参数函数，只需在完整的后验图集上计算它）。使用 pgamma 会给我们一堆落在区间内的概率的后验估计，但我们仍然必须......
...考虑如何将它们聚合成一个估计值（这很容易做到，但对于自动熟悉第一种方法的懒惰分析师来说是额外的思考）。

【解决方案3】：

您可以只使用pgamma 和fit 中的估计参数。

b <- fit$estimate
#   shape     rate 
#1.739679 1.815995 

pgamma(1.5, b[1], b[2]) - pgamma(1, b[1], b[2])
# [1] 0.1896032

谢谢。但是P(x > 2)呢？

查看lower.tail 参数：

pgamma(q, shape, rate = 1, scale = 1/rate, lower.tail = TRUE, log.p = FALSE)

默认情况下，pgamma(q) 计算 Pr(x <= q)。设置lower.tail = FALSE 给出Pr(x > q)。所以你可以这样做：

pgamma(2, b[1], b[2], lower.tail = FALSE)
# [1] 0.08935687

或者你也可以使用

1 - pgamma(2, b[1], b[2])
# [1] 0.08935687

【讨论】：

你太棒了！ lower.tail = TRUE 非常有用。非常感谢。
我不同意这种方法，它忽略了伽马分布参数的最大似然估计周围的后验不确定性。也就是说，如果拟合的参数值完全正确，这种方法会给出正确的答案。但事实并非如此。相反，我们应该将此处计算的概率整合到形状和速率参数的完整联合后验分布上。如果您有兴趣，我可以编写一个快速示例来说明如何执行此操作，但它可能会使用一些不熟悉的包/统计拟合技术。
@ZheyuanLi 是的，你是对的。 TRUE 代表<，这是谈论概率的经典方式。我想知道你是否熟悉主成分分析，我在stackoverflow.com/questions/41022927/… 有关于这个主题的另一个问题。你也能帮我吗？非常感谢。
@ZheyuanLi 抱歉，如果我说你不知道集成是如何工作的。我同意估计的不确定性是已知的，但你和我都同意参数估计中的这种不确定性不会被纳入概率估计Pr(1 < x <= 1.5)，除非我们做一些花哨的事情。我不一定会给自己贴上贝叶斯的标签，但在这种情况下，我喜欢 MCMC 技术带来的不确定性传播的便利性。