机器算法验证 - 加权伯努利随机变量之和的 CDF 是多少？ - 吾爱随笔录

加权伯努利随机变量之和的 CDF 是多少？

机器算法验证分布累积分布函数近似伯努利分布泊松二项分布

2022-03-07 03:37:56

假设我们有一个随机变量定义为个伯努利变量的总和，每个变量都有不同的成功概率和不同的（固定）权重。权重是正数，介于 ~0.001-1,000 之间，有少数异常值。如果它会使问题更易于处理，则可以丢弃非常大和非常小的异常值。 $Y$ $N$ $X_i$ $p_i$ $w_i$

形式上， $Y = \sum X_i w_i$

其中和。 $\Pr(X_i=1)=p_i$ $\Pr(X_i=0)=1-p_i$

$N$ 的范围可以从 ~10-10,000。我想快速计算 $\Pr(Y<=k)$ 的近似值（其中 $k$ 给出）。

我正在处理的数据来自现实生活中的赌注，每个赌注都有隐含赔率 $p_i$ 和赌注 $w_i$ 。以下图表显示了概率和权重的粗略分布：

计算近似值的一种方法是通过蒙特卡罗模拟，但对于较大的 $N$ 值，这可能会变慢。人们也可以使用中心极限定理来计算大（或者如果有少数大权重） $N$ ，准确度很差。最后，有一些方法可以通过使用渐近扩展来改进 CLT 方法，例如： $N$

沃尔科娃，AY（1996 年）。对独立随机指标之和的中心极限定理的改进。概率论及其应用 40, 791-794。

但是，据我所见，精炼的近似值仅指定在所有权重（标准泊松二项分布）。有没有一种方法可以计算出适用于大小的答案，其误差比纯 CLT 方法更小，而无需求助于蒙特卡洛模拟？换句话说，是否可以扩展改进的封闭式近似以允许权重不等于 1？ $Y$ $W_i$ $N$

一些不能完全回答我的相关问题：

CLT 可用于不同伯努利变量的加权和？（纯CLT太不准确了）
如何找到正非整数权重的独立伯努利随机变量的加权和的分布（无答案）
https://math.stackexchange.com/questions/1546366/distribution-of-weighted-sum-of-bernoulli-rvs（[Raghavan , 1988]中提到的边界对于小 N 或计算 CDF 时非常非常宽接近平均值）

4个回答

一种可能性是使用鞍点近似。为此，我们需要 mgf（矩生成函数）及其对数 cgf（累积量生成函数）。带有参数的伯努利变量的 mgf为然后并且 cgf 是那么我们可以按照鞍点逼近如何工作？. $p_i$

M_{i} (t) = E e^{t X_{i}} = (1 - p_{i}) + p_{i} e^{t}

$\DeclareMathOperator{\E}{\mathbb{E}} M_i(t)=\E e^{t X_i}=(1-p_i)+p_i e^t$

E e^{t Y} = \prod_{i} E e^{t w_{i}} = \prod_{i} M_{i} (t w_{i})

$\E e^{tY}=\prod_i \E e^{t w_i}=\prod_i M_i(t w_i)$

K (t) = \sum_{i = 1}^{n} \log (1 - p_{i} + e^{t w_{i}})

$K(t)=\sum_{i=1}^n \log\left( 1-p_i+e^{t w_i} \right)$

使用正态近似

μ = \sum_{i = 1}^{n} p_{i} w_{i}

$\mu = \sum_{i=1}^n p_i w_i$

和方差

σ^{2} = \sum_{i = 1}^{n} w_{i}^{2} p_{i} (1 - p_{i})

$\sigma^2 = \sum_{i=1}^n w_{i}^2 p_i (1-p_i)$

这应该非常准确，尤其是在很大的情况下。 $n$

我的直觉会建议蒙特卡洛的以下变体（尽管我现在想不出证明）。的给定实现开始，计算其并检查是否小于。然后选择一个概率为（只需对一个均匀变量进行采样并使用二等分法在您计算一次的向量中找到它的位置即可轻松高效地完成）开始），从它的分布中重新采样，如果它改变了更新（平均少于一个添加）。再次检查是否小于 $X$ $Y$ $k$ $i$ $w_i/(\Sigma_j w_j)$ cumsum $w_i/(\Sigma_j w_j)$ $X_i$ $Y$ $k$ . 对多次迭代和几个起点（类似于 MCMC 链）重复此方法，以保持小于的次数。该方法类似于 Monte Carlo，但有选择地更频繁地更新元素，这些元素往往对产生更大的影响。 $Y$ $k$ $X_i$ $Y$

尽管您确实要求没有蒙特卡罗模拟的解决方案，但我会强调在这种情况下它的计算成本不会太高，特别是如果您使用并行计算。这是一个R示例

#####
# simulate data
set.seed(91874947)
N  <- 10000 
ps <- ifelse(.4 > runif(N), rbeta(N, 5, 8), rbeta(N, 14, 5)) 
w  <- exp(rnorm(N, -1, .66) + (.2 > runif(N)) * rnorm(N, 1.4, .2))

hist(    ps, breaks = 50)
hist(log(w), breaks = 50)

#####
# define simulation function. `n_sim`` is number of simulations
sim_expr <- function(n_sim = 100000){
  # setup cluster
  require(parallel)
  cl <- makeCluster(7)
  on.exit(stopCluster(cl))
  clusterSetRNGStream(cl)
  clusterExport(cl, c("ps", "w", "N"))

  # run simulation
  parSapply(cl, 1:n_sim, function(...){
    y <- ps > runif(N)
    sum(w * y) 
  })
}

# run simulations and check run time. See the `elapsed` time (time in seconds)
set.seed(37219838)
system.time(s1 <- sim_expr())
#R>   user  system elapsed 
#R>   0.11    0.03   14.22
set.seed(4382482)
system.time(s2 <- sim_expr())
#R>   user  system elapsed 
#R>   0.08    0.00   14.11

# yields almost the same density estimate
plot (density(s1))
lines(density(s2), col = "red")

其它你可能感兴趣的问题

上一篇如何在高斯过程回归中为 RBC 内核选择长度尺度范围？下一篇奥兹国会不会有不开心的部落？