algorithm - Eratosthenes 概率筛法

Question

考虑以下算法。

function Rand():
    return a uniformly random real between 0.0 and 1.0

function Sieve(n):

    assert(n >= 2)

    for i = 2 to n
        X[i] = true

    for i = 2 to n
        if (X[i])
            for j = i+1 to n
                if (Rand() < 1/i)
                    X[j] = false

    return X[n]

Sieve(k) 作为 k 的函数返回 true 的概率是多少？

Answer 1

让我们递归地定义一系列随机变量:

令 X_k,r 表示指标变量，取值 1 当且仅当 X[k] == true变量 i 取值 r 的迭代。

为了使用更少的符号并使代码更直观，我们将只写 X_k,i 这是有效的，尽管自从 i 取值 i 当第一个引用循环中的变量而后者引用变量的值时会造成混淆。

现在我们注意到:

P(X_k,i ~ 0) = P(X_k,i-1 ~ 0) + P(X_k,i-1 ~ 1) * P(X_k-1,i-1 ~ 1) * 1/i

(使用 ~ 代替 = 只是为了使其易于理解，因为 = 否则会有两个不同的含义并且看起来很困惑)。

这个等式成立是因为 X[k] 在 i 迭代结束时为 false 或者因为它在 i-1 结束时为 false，或者在该点为 true，但在最后一次迭代中 X[k-1] 为 true，因此我们进入循环并以 1/i 的概率更改 X[k]。这些事件是互斥的，因此没有交集。

递归的基础就是 P(X_k,1 ~ 1) = 1 和 P(X_2,i ~ 1) = 1 .

最后，我们简单地注意到 P(X[k] == true) = P(X_k,k-1 ~ 1)。

这可以很容易地编程。这是一个使用 memoisation 的 javascript 实现(如果使用嵌套索引比字典索引的字符串连接更好，您可以进行基准测试，您还可以重新设计计算以保持相同的运行时复杂性，但不会通过构建自下而上和不是自上而下)。自然地，该实现将具有 O(k^2) 的运行时复杂度，因此它对于任意大的数字不切实际:

function P(k) {
   if (k<2 || k!==Math.round(k)) return -1;
   var _ = {};
   function _P(n,i) {
      if(n===2) return 1;
      if(i===1) return 1;
      var $ = n+'_'+i;
      if($ in _) return _[$];
      return _[$] =  1-(1-_P(n,i-1) + _P(n,i-1)*_P(n-1,i-1)*1/i);
   }
   return _P(k,k-1);
}

P(1000); // 0.12274162882390949

更有趣的是 1/i 概率如何改变事物。 IE。概率是否收敛到 0 或某个其他值，如果收敛，改变 1/i 会如何影响它。

当然，如果您在 mathSE 上提问，您可能会得到更好的答案 - 这个答案非常简单，我相信有一种方法可以对其进行操作以获得直接公式。