haskell - Haskell 中的随机数采样序列

Question

我需要小列表的高斯随机数进行模拟，因此我尝试了以下操作:

import System.Random

seed = 10101
gen = mkStdGen seed

boxMuller mu sigma (r1,r2) =  mu + sigma * sqrt (-2 * log r1) * cos (2 * pi * r2) 

这只是 Box-Muller 算法 - 给定 [0,1] 区间内的 r1, r2 均匀随机数，它返回一个高斯随机数。

normals 0 g = [] 
normals n g = take n $ map (boxMuller 0 1) $ pairs $ randoms g
    where pairs (x:y:zs) = (x,y):(pairs zs)

因此，每当我需要随机数列表时，我都会使用这个 normals 函数。

问题一定很明显:它总是生成相同的序列，因为我总是使用相同的种子!我没有得到新的序列，我一直只得到序列的前 n 个值。

当我输入时，我明显假装的是:

x = normal 10 
y = normal 50

我将 x 作为 map (boxMuller 0 1) $pairs $randoms g 的前 10 个值，将 y 作为此列表中接下来的 50 个值，依此类推。

当然这是不可能的，因为给定相同的输入，函数必须始终返回相同的值。如何摆脱这个陷阱？

Answer 1

我认为，在抽象生成器的 monad 中进行需要随​​机数的计算是最干净的事情。该代码如下所示:

我们将把 StdGen 实例放入状态 monad 中，然后在状态 monad 的 get 和 set 方法上提供一些糖，以给我们随机数。

首先，加载模块:

import Control.Monad.State (State, evalState, get, put)
import System.Random (StdGen, mkStdGen, random)
import Control.Applicative ((<$>))

(通常我可能不会指定导入，但这可以很容易地理解每个函数的来源。)

然后我们将定义“需要随机数的计算”monad；在本例中，State StdGen 的别名称为 R。 (因为“随机”和“兰德”已经意味着其他含义。)

type R a = State StdGen a

R 的工作方式是:定义一个需要随机数流(单子(monad)“ Action ”)的计算，然后使用 runRandom“运行它”:

runRandom :: R a -> Int -> a
runRandom action seed = evalState action $ mkStdGen seed

这需要一个操作和一个种子，并返回操作的结果。就像通常的evalState、runReader等

现在我们只需要 State monad 周围的糖。我们使用get来获取StdGen，并使用put来安装新状态。这意味着，要获得一个随机数，我们可以这样写:

rand :: R Double
rand = do
  gen <- get
  let (r, gen') = random gen
  put gen'
  return r

我们获取随机数生成器的当前状态，用它来获取新的随机数和新的生成器，保存随机数，安装新的生成器状态，然后返回随机数。

这是一个可以使用 runRandom 运行的“操作”，所以让我们尝试一下:

ghci> runRandom rand 42
0.11040701265689151                           
it :: Double     

这是一个纯函数，因此如果您使用相同的参数再次运行它，您将得到相同的结果。杂质保留在您传递给 runRandom 的“ Action ”内。

无论如何，您的函数需要随机数对，所以让我们编写一个操作来生成一对随机数:

randPair :: R (Double, Double)
randPair = do
  x <- rand
  y <- rand
  return (x,y)

使用 runRandom 运行此命令，您将看到这对数字中的两个不同数字，正如您所期望的那样。但请注意，您不必为“rand”提供参数；这是因为函数是纯函数，但“rand”是一个 Action ，不需要是纯函数。

现在我们可以实现您的法线函数。 boxMuller 正如您在上面定义的那样，我只是添加了一个类型签名，这样我就可以了解正在发生的事情，而无需阅读整个函数:

boxMuller :: Double -> Double -> (Double, Double) -> Double
boxMuller mu sigma (r1,r2) =  mu + sigma * sqrt (-2 * log r1) * cos (2 * pi * r2)

实现所有辅助函数/操作后，我们终于可以编写 normals，这是一个 0 个参数的操作，返回一个(延迟生成的)正态分布 double 的无限列表:

normals :: R [Double]
normals = mapM (\_ -> boxMuller 0 1 <$> randPair) $ repeat ()

如果您愿意，您也可以写得更简洁:

oneNormal :: R Double
oneNormal = do
    pair <- randPair
    return $ boxMuller 0 1 pair

normals :: R [Double]
normals = mapM (\_ -> oneNormal) $ repeat ()

repeat () 为单子(monad) Action 提供了一个无限的流来调用函数(这也是法线结果无限长的原因)。我最初编写了[1..]，但我重写了它，从程序文本中删除了无意义的1。我们不是对整数进行操作，我们只是想要一个无限列表。

总之，就这样吧。要在实际程序中使用它，您只需在 R Action 中执行需要随机法线的工作即可:

someNormals :: Int -> R [Double]
someNormals x = liftM (take x) normals

myAlgorithm :: R [Bool]
myAlgorithm = do
   xs <- someNormals 10
   ys <- someNormals 10
   let xys = zip xs ys
   return $ uncurry (<) <$> xys

runRandom myAlgorithm 42

应用一元 Action 编程的常用技术；在 R 中保留尽可能少的应用程序，这样事情就不会太困惑。

哦，还有另一件事:惰性可以干净地“泄漏”到 monad 边界之外。这意味着你可以写:

take 10 $ runRandom normals 42

它会起作用。