随机算法表现不如预期

Question

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

使用 log（log n）位维护计数器X

• 将计数器X初始化为0

• 当一个项目到达时，以概率（½）^X将X增加1

• 当流结束时，输出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]

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

我在 Matlab 中测试了相同的算法，它在那里运行良好，所以我认为要么是：

1.我的随机数生成器存在问题；或者 2.在Double中将1/2提高到一个大的n

请建议下一步操作。

- xiaolingxiao

---

如果在Matlab中执行正常，很可能不是算法的问题。我不知道这是什么语言，但你应该在那个stackoverflow上发布问题。 - ElKamina

你的数学符号在我的手机上显示不清楚。你能使用ASCII码，或者至少更常见的符号吗？ - dfeuer

实际上 0.5^2000 :: Double 是零，但我看不出这会在这里引起麻烦。 - chi

这不是你的问题，但请注意像这样传递StdGen是容易出错的，因为很容易使用旧的或两次使用新的。话虽如此，就我所见，您的代码似乎正确地传递了它们。为了防止这些陷阱，在将来考虑使用Control.Monad.Random中的Rand单子。 - chi

dfeuer我刚刚做了。如果我一直在传递stdGen，使用Rand会有什么区别吗？ - xiaolingxiao

Answer 1

问题其实很简单：使用 randomR（1,100） 会排除第一个百分比范围内的值，因此在 1/2 的高次幂处完全被截断（这些值都位于该小区间内）。实际上是一个普遍问题：范围应该从零开始，而不是从一开始^†，除非有特殊原因。

但是，为什么要在第一个地方使用100个范围呢？我只需让它变成

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

---

^†_{我知道，Matlab在很多地方都会出错。这只是关于该语言的许多可怕之处之一。（参考资料）}

---

与您的问题无关：正如 chi 所指出的那样，如果您使用适当的随机单子（random monad），而不是手动传递 StdGen，此类代码看起来会好看得多。

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

- leftaroundabout

---

现在有没有一种方法可以概括 tos 的签名，即 MonadState Int m => m Int vs State Int Int？我似乎找不到一些 MonadRVar。 - xiaolingxiao

@chibro2：你可以随时使用lift来访问在RVal之上的变换器堆栈，或者你可以在堆栈顶部直接使用RValT。 - leftaroundabout