在Haskell中随机数的抽样序列 [英] Sampling sequences of random numbers in Haskell

问题描述

  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算法 - 给定r1，r2 uniform uniform

  normals 0 g = [] 
 [0,1]间隔中的数字将返回一个高斯随机数。法线ng =取n $ map（boxMuller 0 1）$ pairs $ randoms g 
其中pairs（x：y：zs）=（x，y）:(对zs）
  normals 函数。  
 
 
这个问题一定是显而易见的：它总是产生相同的序列，因为我总是使用相同的种子！我没有得到新的序列，我只是一直得到序列的前n个值。
 
 
 
我当时假装清楚的是，当我输入：
  x =正常10 
y =正常50 
  
我有x是 map的前10个值（boxMuller 0 1）$ pairs $ randoms g 和y作为该列表中的下50个值，依此类推。  
 
 
当然这是不可能的，因为一个函数在给定相同的输入时总是返回相同的值。我怎么逃避这个陷阱？
解决方案
我认为，你的计算需要在monad中抽取随机数是最干净的事情。以下是代码的样子：
 
 
 
我们将StdGen实例置于状态monad中，然后提供状态monad的get和set方法的一些参数首先，加载模块： 
 
 
  导入Control.Monad.State（State，evalState，get，put）
导入System.Random（StdGen，mkStdGen，随机）
导入Control.Applicative（（< $>））
  
（通常我可能不会指定导入，但这可以很容易地理解每个函数的来源。 ）
 
 
然后我们将定义我们的需要随机数的计算monad;在这种情况下， State StdGen 的别名称为 R 。 （因为随机和兰德已经意味着别的东西了。）
  type R a =状态StdGen a 
  
 R的工作方式是：定义一个计算需要一串随机数（monadicaction），然后你用 runRandom 运行它：
  runRandom :: R a  - > Int  - > a 
 runRandom action seed = evalState action $ mkStdGen seed 
  
这需要一个动作和一个种子并返回操作的结果。就像通常的 evalState ， runReader 等。 
 
 
现在我们只需要在州monad周围加糖。我们使用 get 来获取StdGen，并且使用 put 来安装一个新状态。这意味着，要得到一个随机数，我们会写：$​​ b 
 $ b 
  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是一个动作，它不一定是纯粹的。
 
 
 
现在我们可以实现 normals 函数。  boxMuller 就像你上面定义的那样，我只是添加了一个类型签名，这样我就可以理解正在发生的事情，而无需阅读整个函数： 
 
 
  boxMuller :: Double  - >双 - > （Double，Double） - > Double 
 boxMuller mu sigma（r1，r2）= mu + sigma * sqrt（-2 * log r1）* cos（2 * pi * r2）
  
执行完所有帮助程序函数/操作后，我们最终可以编写 normals ，一个0参数的操作，返回一个（lazily-generated）正常分布的双精度型无限列表： 
 
 
  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（）给出了monadic行动无限的无用的任何东西来调用函数（并且是法线无限长的结果）。我最初编写了 [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 
  
，它会起作用。
 
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) 
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)
So I used this normals function everytime I needed my list of random numbers. 

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.

What I was pretending clearly was that, when I type:
x = normal 10 
y = normal 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?
解决方案
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:

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.

First, load the modules:
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.)

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
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
This takes an action and a seed, and returns the results of the action.  Just like the usual evalState, runReader, etc. 

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.

This is an "action" that can be run with runRandom, so let's try it:
ghci> runRandom rand 42
0.11040701265689151                           
it :: Double     
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)
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.

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)
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 () 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.

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

Oh, and on other thing: laziness can "leak" outside of the monad boundary cleanly.  This means you can write:
take 10 $ runRandom normals 42
and it will work.

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