伯努利分布式布尔值的 Haskell 无限列表

Question

我需要一个有偏见的随机布尔值列表。每个布尔值都需要具有相同的真概率（伯努利分布）。这些布尔值被传递给一个函数，该函数为每个输入布尔值生成零个或多个输出布尔值。我需要一个无限列表，因为我事先不知道需要多少布尔值才能提供足够的输出。请参阅以下（简化）代码：

import System.Random.MWC
import System.Random.MWC.Distributions

foo :: [Bool] -> [Bool] -- foo outputs zero or more Bools per input Bool

main = do
  gen <- create
  bits <- sequence . repeat $ bernoulli 0.25 gen
  print . take 32 . foo $ bits

不幸的是，这段代码只挂在main. 我猜在某处发生了一些非懒惰的事情Control.Monad.ST？

（我可以用 做这样的事情System.Random.randoms，但结果值没有所需的分布。）

我可以在继续使用System.Random.MWC图书馆的同时解决这个问题吗？或者这是否需要我切换到替代实现？

Answer 1

该mwc-random包提供了两个PrimMonad实例，一个 forIO和另一个 for ST s。只要ST计算在所有状态标签上参数化s，我们就可以运行计算并使用 提取值runST :: (forall s. ST s a) -> a。这本身并不是很有用，因为我们会丢失状态：随机生成器的种子，但 mwc-random 还提供了处理种子的明确方法：

save :: PrimMonad m => Gen (PrimState m) -> m Seed
restore :: PrimMonad m => Seed -> m (Gen (PrimState m))

只要生成器在forall s. ST s.

{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}

import System.Random.MWC
import Control.Monad.ST
import System.Random.MWC.Distributions

randomStream :: forall s a. (forall s. GenST s -> ST s a) -> GenST s -> ST s [a]
randomStream item = go
    where
        go :: forall s. GenST s -> ST s [a]
        go gen = do
            x <- item gen
            seed <- save gen
            return (x:runST (restore seed >>= go))

有了这个，我们可以把你的例子写成

main = do
    bits <- withSystemRandom (randomStream (bernoulli 0.25))
    print . take 32 $ bits

我们实际上可以构建比对流中的每个项目使用相同的生成器更复杂的生成器。我们可以沿着流线程化一个状态，这样每个值都可以依赖于前一个值的结果。

unfoldStream :: forall s a b. (forall s. b -> GenST s -> ST s (a, b)) -> b -> GenST s -> ST s [a]
unfoldStream item = go
    where
        go :: forall s. b -> GenST s -> ST s [a]
        go b gen = do
            (x,b') <- item b gen
            seed <- save gen
            return (x:runST (restore seed >>= go b'))

以下示例流的结果在每次结果为 时增加可能性False。

import Control.Monad.Primitive

interesting :: (PrimMonad m) => Double -> Gen (PrimState m) -> m (Bool, Double)
interesting p gen = do
    result <- bernoulli p gen
    let p' = if result then p else p + (1-p)*0.25
    return (result, p')

main = do
    bits <- withSystemRandom (unfoldStream interesting 0)
    print . take 32 $ bits