随机算法未按预期运行

Posted 2023-03-06

【中文标题】随机算法未按预期运行

【英文标题】：Randomized algorithm not behaving as expected

【发布时间】：2015-12-06 04:49:44

【问题描述】：

我正在实现一个近似计数算法，我们：

使用 log (log n) bits 维护一个计数器 X

初始化X为0

当物品到达时，以 (½) 的概率将 X 增加 1^X

当流结束时，输出 2^X - 1 使得 E[2^X]= n + 1

我的实现如下：

import System.Random

type Prob   = Double
type Tosses = Int

-- * for sake of simplicity we assume 0 <= p <= 1
tos :: Prob -> StdGen -> (Bool,StdGen)
tos p s = (q <= 100*p, s')
  where (q,s') = randomR (1,100) s

toses :: Prob -> Tosses -> StdGen -> [(Bool,StdGen)]
toses _ 0 _ = []
toses p n s = let t@(b,s') = tos p s in t : toses p (pred n) s'

toses' :: Prob -> Tosses -> StdGen -> [Bool]
toses' p n = fmap fst . toses p n

morris :: StdGen -> [a] -> Int
morris s xs = go s xs 0 where
  go _ []     n = n
  go s (_:xs) n = go s' xs n' where
    (h,s') = tos (0.5^n) s 
    n'     = if h then succ n else n

main :: IO Int
main = do
  s <- newStdGen
  return $ morris s [1..10000]

问题是我的 X 对于任何 |stream| > 2 总是不正确，而且似乎对于所有 StdGen 和 |stream| > 1000、X = 7 都是不正确的

我在 Matlab 中测试了相同的算法，它在那里工作，所以我假设它要么是

我的随机数生成器有问题，或者

在Double中将 1/2 提高到一个大的 n

请提出前进的道路？

【问题讨论】：

如果它在 Matlab 中工作，它不可能是算法问题。我不知道这是什么语言，但你应该在那个 *** 中发帖

你的数学符号在我的手机上显示不好。你能使用 ASCII，或者至少是更常用的符号吗？

确实0.5^2000 :: Double 为零，但我看不出这会不会在这里造成麻烦。

这不是您的问题，但请注意，像这样传递StdGen 很容易出错，因为使用旧的或使用新的很容易两次。话虽如此，据我所知，您的代码似乎正确地传递了它们。为防止此类陷阱，将来，请考虑使用 Rand 之类的 monad，来自 Control.Monad.Random。

我刚刚做了。和 @chi 如果我传递 stdGen 会使用 Rand 有所作为吗？

【参考方案1】：

问题实际上非常简单：使用randomR (1,100)，您排除了第一个百分比内的值，因此您在 1/2 的高次幂处完全截止（所有这些都位于那个小区间内）。实际上是一个普遍的事情：ranges should start at zero，不是一个^†，除非有特定的原因。

但是为什么一开始就使用 100 的范围呢？我会成功的

tos :: Prob -> StdGen -> (Bool,StdGen)
tos p s = (q <= p, s')
  where (q,s') = randomR (0,1) s

---

^†_{我知道，Matlab 到处都是这个错误。只是many 那门语言的可怕之处之一。}

---

与您的问题无关：正如 chi 所说，如果您使用合适的随机 monad，这种代码看起来会更好，而不是手动传递 StdGens。

import Data.Random
import Data.Random.Source.Std

type Prob   = Double

tos :: Prob -> RVar Bool
tos p = do
  q <- uniform 0 1
  return $ q <= p

morris :: [a] -> RVar Int
morris xs = go xs 0 where
  go []     n = return n
  go (_:xs) n = do
    h <- tos (0.5^n)
    go xs $ if h then succ n else n

morrisTest :: Int -> IO Int
morrisTest n = do
  runRVar (morris [1..n]) StdRandom

【讨论】：

现在有没有办法泛化tos的签名，即MonadState Int m => m Int vs State Int Int？我似乎找不到一些MonadRVar

@chibro2：您始终可以使用lift 进入RVal 顶部的变压器堆栈，或者您可以在堆栈顶部使用RValT 本身。