18

我需要一个小的高斯随机数列表来进行模拟，所以我尝试了以下方法：

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 成为该列表的前 10 个值，map (boxMuller 0 1) $ pairs $ randoms gy 成为该列表中的下 50 个值，依此类推。

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

4

3 回答 3

28

我认为在抽象生成器的单子中进行需要随​​机数的计算是最干净的事情。这是该代码的样子：

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

首先，加载模块：

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

（通常我可能不会指定导入，但这很容易理解每​​个函数的来源。）

然后我们将定义我们的“需要随机数的计算”单子；在这种情况下，为State StdGencalled的别名R。（因为“随机”和“兰德”已经有别的意思了。）

type R a = State StdGen a

R 的工作方式是：定义一个需要随机数流的计算（一元“动作”），然后使用以下命令“运行它” 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 的“动作”中。

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

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

用 runRandom 运行它，你会看到这对中的两个不同的数字，正如你所期望的那样。但请注意，您不必提供带有参数的“rand”；那是因为函数是纯的，但“rand”是一个动作，它不必是纯的。

现在我们可以实现你的normals功能了。 boxMuller正如你在上面定义的那样，我只是添加了一个类型签名，这样我就可以理解发生了什么而不必阅读整个函数：

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

实现了所有辅助函数/动作后，我们终于可以编写normals0 个参数的动作，它返回一个（延迟生成的）正态分布双精度数的无限列表：

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 ()给一元动作一个无限的无源流来调用函数（并且是使法线的结果无限长的原因）。我最初是写[1..]的，但我重写了它以1从程序文本中删除无意义的内容。我们不是对整数进行操作，我们只想要一个无限列表。

无论如何，就是这样。要在实际程序中使用它，您只需在 R 动作中完成需要随机法线的工作：

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

适用于编程一元动作的常用技术；尽可能少地保留您的应用程序R，并且事情不会太混乱。

哦，还有另一件事：懒惰可以干净地“泄漏”到单子边界之外。这意味着您可以编写：

take 10 $ runRandom normals 42

它会起作用。

于 2010-01-24T00:16:55.367 回答
6

您从中获得的列表randoms是无限的，当您使用有限前缀时，您无需丢弃无限尾。您可以将随机数作为附加参数传入，并将未使用的随机数作为附加结果返回，或者您可以将无限的随机数序列存放在 state monad 中。

对于需要提供唯一符号的编译器和其他代码，也会出现类似的问题。这在 Haskell 中只是一个真正的痛苦，因为您在整个代码中线程化状态（随机数生成器或符号生成器）。

我已经使用显式参数和单子完成了随机算法，但没有一个是真正令人满意的。如果你 grok monads，我可能会稍微建议使用一个包含尚未使用的随机数的无限列表的状态 monad。

于 2010-01-21T16:10:28.643 回答
1

你也可以通过使用来回避这个问题newStdGen，然后你每次都会得到一个不同的种子（实际上）。

于 2010-01-21T19:00:06.100 回答