用 R 中的伽马分布估计均值和置信区间的标准差答案

【问题标题】：Estimating the standard deviation from mean and confidence intervals with a gamma distribution in R用 R 中的伽马分布估计均值和置信区间的标准差
【发布时间】：2020-05-10 08:00:38
【问题描述】：

我想在 R 中解决以下问题并将其应用于更大的工作流程。我需要从已知均值和 95% 置信区间的 gamma 分布中估计标准差。

state = c("group1", "group2", "group3")
mean = c(0.970, 0.694, 0.988)
lowers = c(0.527, 0.381, 0.536)
uppers = c(1.87, 1.37, 1.90)

df = data.frame(state=state, mean=mean, lower=lower, upper=upper)

使用 excel 和“求解器”工具，我可以调整标准偏差，以最小化目标 2.5（下限）和 97.5（上限）分布与实际值之间的平方差之和。挑战在于这需要扩展到相当大的数据集并在我的 R 数据框工作流程中进行操作。任何想法如何解决这个问题？

【问题讨论】：

这三个样本的“真相”你知道吗？
我很好奇：如果您有大量数据，那么想出 95% 置信区间比标准差更容易吗？方差可以加法计算，因此对此的分布式计算相对简单。分位数很难在规模上完美地完成，因为它们通常需要排序的数据。
我有大量数据需要应用这个问题。我无权访问用于计算平均值和 95% CI 的基础数据。
好的，有道理。到目前为止，您能够完成什么，即使您认为它效率低下？我的猜测是它可能使用optim 或optimize。
R 中不多。我有一个 Excel 工作簿可以做到这一点 - 一个一个 - 这是可管理的，但肯定不理想。真的希望得到一个代码解决方案，它可以成为我们常规工作流程的一部分并且可重现。我开始看 optim，但很快就迷路了。

标签： r

【解决方案1】：

我认为这个问题最终是一个优化问题，一次处理一行数据。不过，由于您想对其进行缩放，这里有一个近似值，用于查找分布核心参数。

此过程不是优化：它扩展了可能的k（形状）参数的定义范围，并找到最接近您的鞋面的形状/比例组合（给定您的平均值）和较低的分位数。您可以控制k 的粒度，这与您获得容差（这将适合优化）一样好。

因此，这个过程将是不完美的。我希望它能为您提供一个足够快的过程，以获得足够好的准确性。

我将首先演示一个每次运行一行的过程，如guesser1。然后我将扩展它以对任意数量的mean、lower 和upper 执行相同的操作。

已知答案的数据

但首先，我想生成我自己的样本，以便我们知道这个猜测者的“真相”。

library(dplyr)
set.seed(42)
n <- 4
randks <- tibble(
  k = sample(1:10, size = n, replace = TRUE),
  scale = sample(seq(0.5, 2.5, by = 0.5), size = n, replace = TRUE)
) %>%
  mutate(
    samp = map2(k, scale, ~ rgamma(1000, shape = .x, scale = .y)),
    theor_mu = k*scale,
    mu = map_dbl(samp, ~ mean(.x)),
    lwr = map_dbl(samp, ~ quantile(.x, 0.025)),
    upr = map_dbl(samp, ~ quantile(.x, 0.975))
  ) %>%
  select(-samp)
randks
# # A tibble: 4 x 6
#       k scale theor_mu    mu   lwr   upr
#   <int> <dbl>    <dbl> <dbl> <dbl> <dbl>
# 1    10   2       20   19.9   9.47 33.7 
# 2    10   1.5     15   15.1   7.36 25.0 
# 3     3   2        6    5.85  1.08 14.5 
# 4     9   0.5      4.5  4.51  1.99  7.72

Guesser1

一次一行：

guesser1 <- function(mu, lwr, upr, k.max = 10, k.by = 0.01) {
  stopifnot(length(mu) == 1, length(lwr) == 1, length(upr) == 1)
  ks <- seq(0, k.max, by = k.by)[-1]
  L <- sapply(ks, function(k) qgamma(0.025, shape = k, scale = mu / k))
  U <- sapply(ks, function(k) qgamma(0.975, shape = k, scale = mu / k))
  dists <- sqrt((L-lwr)^2 + (U-upr)^2)
  ind <- which.min(dists)
  data.frame(
    k     = ks[ind],
    scale = mu / ks[ind],
    dist  = min(dists),
    lwr   = L[ind],
    upr   = U[ind]
  )
}

在行动：

out1 <- do.call(rbind, Map(guesser1, randks$mu, randks$lwr, randks$upr))
cbind(subset(randks, select = -theor_mu), out1)
#    k scale    mu  lwr   upr     k scale  dist  lwr   upr
# 1 10   2.0 19.88 9.47 33.67 10.00 1.988 0.304 9.53 33.97
# 2 10   1.5 15.06 7.36 25.02 10.00 1.506 0.727 7.22 25.73
# 3  3   2.0  5.85 1.08 14.50  2.76 2.120 0.020 1.10 14.50
# 4  9   0.5  4.51 1.99  7.72  9.55 0.472 0.142 2.12  7.79
###  \____ randks __________/  \____ guessed ____________/

肯定存在一些差异，强调了我最初的断言，即这是不完美的。

猜测者

一次所有行。这是函数中的更多工作，因为它处理矩阵而不仅仅是向量。没问题，我只是想在尽情享受之前一次证明一下。

guessers <- function(mu, lwr, upr, k.max = 10, k.by = 0.01, include.size = FALSE) {
  stopifnot(length(mu) == length(lwr), length(mu) == length(upr))
  # count <- length(mu)
  ks <- seq(0, k.max, by = k.by)[-1]
  # 'ks'      dims: [mu]
  L <- outer(mu, ks, function(m, k) qgamma(0.025, shape = k, scale = m / k))
  U <- outer(mu, ks, function(m, k) qgamma(0.975, shape = k, scale = m / k))
  # 'L' & 'U' dims: [mu, ks]
  dists <- sqrt((L - lwr)^2 + (U - upr)^2)
  inds <- apply(dists, 1, which.min)
  mindists <- apply(dists, 1, min)
  i <- seq_along(mu)
  out <- data.frame(
    k = ks[inds],
    scale = mu / ks[inds],
    dist = mindists,
    lwr = L[cbind(i, inds)],
    upr = U[cbind(i, inds)]
  )
  size <- if (include.size) {
    message("guessers memory: ",
            object.size(list(ks, L, U, dists, inds, mindists, i, out)))
  }
  out
}

在行动：

outs <- guessers(randks$mu, randks$lwr, randks$upr, include.size = TRUE)
# guessers memory: 106400
cbind(subset(randks, select = -theor_mu), outs)
#    k scale    mu  lwr   upr     k scale  dist  lwr   upr
# 1 10   2.0 19.88 9.47 33.67 10.00 1.988 0.304 9.53 33.97
# 2 10   1.5 15.06 7.36 25.02 10.00 1.506 0.727 7.22 25.73
# 3  3   2.0  5.85 1.08 14.50  2.76 2.120 0.020 1.10 14.50
# 4  9   0.5  4.51 1.99  7.72  9.55 0.472 0.142 2.12  7.79
###  \____ randks __________/  \____ guessed (same) _____/

（我在其中包含了一条内存消息，只是为了跟踪它可以扩展多少。现在还不错，而且这个参数绝对不应该在生产中使用。仅供参考。）

比较

使用microbenchmark，我们将每个操作重复几次并比较它们的运行时间。

microbenchmark::microbenchmark(
  g1 = Map(guesser1, randks$mu, randks$lwr, randks$upr),
  gs = guessers(randks$mu, randks$lwr, randks$upr)
)
# Unit: milliseconds
#  expr  min   lq mean median   uq   max neval
#    g1 27.3 28.8 33.9   29.7 33.0 131.1   100
#    gs 13.3 13.6 14.4   13.8 14.6  20.3   100

毫不奇怪，一次性guessers 更快一点。什么时候不会出现这种情况？当行数变得如此之大以至于内存消耗成为问题时。我不知道那是什么。

让我们尝试同样的事情，但将 randks 从 4 行增加到 1000 行并重复基准测试。

n <- 1000
# randks <- ...
nrow(randks)
# [1] 1000
microbenchmark::microbenchmark(
  g1 = Map(guesser1, randks$mu, randks$lwr, randks$upr),
  gs = guessers(randks$mu, randks$lwr, randks$upr),
  times = 10
)
# Unit: seconds
#  expr  min   lq mean median   uq   max neval
#    g1 8.50 8.99 9.59   9.31 9.64 11.72    10
#    gs 3.35 3.44 3.61   3.63 3.77  3.93    10

所以它肯定更快。 1000 次估计的中位运行时间为 3.63 秒，因此它似乎完成了大约 300/秒。

guessers memory: 24066176

(24 MiB) 实际上，这似乎一点也不坏。减少 k.by 以提高准确性，或增加 k.by 以加快速度。

【讨论】：