伯努利分布布尔值的 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 实例，一个用于IO，另一个用于ST s。只要在所有状态标签s 上参数化ST 计算，我们就可以运行计算并使用runST :: (forall s. ST s a) -> a 提取值。这本身不会很有用，因为我们会丢失状态：随机生成器的种子，但 mwc-random 也提供了handle the seeds 的显式方法：

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

Answer 2

罪魁祸首是sequence . repeat - 这将挂起（几乎？）每个 monad，因为您必须执行可能无限数量的效果。

最简单的解决方案是使用不同的库 - 如果您依赖 mwc-random 生成的数字的质量，这可能是不可能的。下一个最简单的解决方案是将foo 重写为具有[IO Bool] -> IO [Bool] 类型并将其传递给repeat (bernoulli 0.25 gen) - 这将允许foo 选择何时停止执行由无限列表产生的效果。但是在IO 中包含你的逻辑并不是很好。

当您需要无限的随机数列表时，标准技巧是使用纯函数f :: StdGen -> (Result, StdGen)。然后unfoldr (Just . f) :: StdGen -> [Result]，输出是一个无限列表。乍一看，mwc-random 似乎只有一元函数，没有纯接口。然而，事实并非如此，因为ST s 是PrimMonad 的一个实例。您还拥有将Gen 转换为Seed 的函数。使用这些，您可以获得任何一元的纯 RNG 函数：

{-# LANGUAGE RankNTypes #-}

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

pureRand :: (forall s . GenST s -> ST s t) -> Seed -> (t, Seed) 
pureRand f s = runST $ do 
  s'  <- restore s
  r   <- f s' 
  s'' <- save s' 
  return (r, s'')

pureBernoulli :: Double -> Seed -> (Bool, Seed)
pureBernoulli a = pureRand (bernoulli a) 

foo :: [Bool] -> [Bool]
foo = id 

main = do
  gen <- create >>= save
  let bits = unfoldr (Just . pureBernoulli 0.25) gen 
  print . take 32 . foo $ bits

不幸的是，mwc-random 默认不公开这种接口，但很容易获得。

另一种选择更可怕——使用不安全的函数。

import System.IO.Unsafe

repeatM rand = go where
  go = do
    x  <- rand
    xs <- unsafeInterleaveIO go
    return (x : xs)

main2 = do
  gen <- create
  bits <- repeatM (bernoulli 0.25 gen) 
  print . take 32 . foo $ bits

当然，这伴随着unsafe 周围的常见警告 - 仅当您对纯函数非常不便时才使用它。 unsafeInterleaveIO 可能会重新排序或从不执行效果 - 例如，如果 foo 忽略一个元素，它将永远不会被计算并且更新存储在 gen 中的状态的相应效果可能不会发生。例如，以下内容将不打印任何内容：

snd <$> ((,) <$> unsafeInterleaveIO (putStrLn "Hello") <*> return ())