Haskell中随机数的采样序列

问题描述

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

I need small lists of gaussian random numbers for a simulation and so I tried the following:

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 均匀随机数，它返回一个高斯随机数.

This is just the Box-Muller algorithm - given r1, r2 uniform random numbers in the [0,1] interval it returns a gaussian random number.

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 函数.

So I used this normals function everytime I needed my list of random numbers.

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

The problem with that must be apparent: it generates always the same sequence cause I'm using always the same seed! I'm not getting new sequences, I'm only getting the first n values of the sequence all the time.

当我输入:

x = normal 10
y = normal 50

我希望 x 是 map (boxMuller 0 1) $pairs $ randoms g 的前 10 个值，y 是此列表中的下一个 50 个值，依此类推.

I'd have x to be the first 10 values of map (boxMuller 0 1) $ pairs $ randoms g and y to be the next 50 values on this list, and so on.

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

Of course this is impossible, cause a function must always return the same values given the same input. How do I escape this trap?

推荐答案

我认为在抽象掉生成器的 monad 中进行需要随​​机数的计算是最干净的事情.下面是代码的样子:

I think that doing your computations that require random numbers inside a monad that abstracts away the generator is the cleanest thing. Here is what that code would look like:

我们将把 StdGen 实例放在一个 state monad 中，然后在 state monad 的 get 和 set 方法上提供一些糖来为我们提供随机数.

We are going to put the StdGen instance in a state monad, then provide some sugar over the state monad's get and set method to give us random numbers.

首先，加载模块:

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

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

(Normally I would probably not specify the imports, but this makes it easy to understand where each function is coming from.)

然后我们将定义我们的需要随机数的计算"monad；在这种情况下，State StdGen 的别名称为 R.(因为随机"和兰德"已经意味着别的东西了.)

Then we will define our "computations requiring random numbers" monad; in this case, an alias for State StdGen called R. (Because "Random" and "Rand" already mean something else.)

type R a = State StdGen a

R 的工作方式是:定义一个需要随机数流的计算(一元动作")，然后使用 runRandom 来运行"它:

The way R works is: you define a computation that requires a stream of random numbers (the monadic "action"), and then you "run it" with runRandom:

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

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

This takes an action and a seed, and returns the results of the action. Just like the usual evalState, runReader, etc.

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

Now we just need sugar around the State monad. We use get to get the StdGen, and we use put to install a new state. This means, to get one random number, we would write:

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

我们获取随机数生成器的当前状态，用它来获取一个新的随机数和一个新的生成器，保存随机数，安装新的生成器状态，并返回随机数.

We get the current state of the random number generator, use it to get a new random number and a new generator, save the random number, install the new generator state, and return the random number.

这是一个可以用runRandom运行的动作"，让我们试试看:

This is an "action" that can be run with runRandom, so let's try it:

ghci> runRandom rand 42
0.11040701265689151
it :: Double

这是一个纯函数，所以如果你用相同的参数再次运行它，你会得到相同的结果.杂质留在您传递给 runRandom 的动作"中.

This is a pure function, so if you run it again with the same arguments, you will get the same result. The impurity stays inside the "action" that you pass to runRandom.

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

Anyway, your function wants pairs of random numbers, so let's write an action to produce a pair of random numbers:

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

用 runRandom 运行这个，你会看到这对中的两个不同的数字，正如你所期望的.但是请注意，您不必为rand"提供参数；那是因为函数是纯的，而rand"是一个动作，不需要是纯的.

Run this with runRandom, and you'll see two different numbers in the pair, as you'd expect. But notice that you didn't have to supply "rand" with an argument; that's because functions are pure, but "rand" is an action, which need not be pure.

现在我们可以实现您的 normals 功能.boxMuller 就像你在上面定义的那样，我只是添加了一个类型签名，所以我无需阅读整个函数就可以理解发生了什么:

Now we can implement your normals function. boxMuller is as you defined it above, I just added a type signature so I can understand what's going on without having to read the whole function:

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

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

With all the helper functions/actions implemented, we can finally write normals, an action of 0 arguments that returns a (lazily-generated) infinite list of normally-distributed Doubles:

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

如果你愿意，你也可以不那么简洁地写:

You could also write this less concisely if you want:

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

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

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

repeat () gives the monadic action an infinite stream of nothing to invoke the function with (and is what makes the result of normals infinitely long). I had initially written [1..], but I rewrote it to remove the meaningless 1 from the program text. We are not operating on integers, we just want an infinite list.

无论如何，就是这样.要在实际程序中使用它，您只需在 R 操作中执行需要随机法线的工作:

Anyway, that's it. To use this in a real program, you just do your work requiring random normals inside an 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

适用于编程一元动作的常用技术；尽可能少地在 R 中保留您的应用程序，这样事情就不会太乱了.

The usual techniques for programming monadic actions apply; keep as little of your application in R as possible, and things won't be too messy.

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

Oh, and on other thing: laziness can "leak" outside of the monad boundary cleanly. This means you can write:

take 10 $ runRandom normals 42

它会起作用.

这篇关于Haskell中随机数的采样序列的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！