secuencias de muestreo de números aleatorios en Haskell

Question

que necesitan pequeñas listas de números aleatorios gaussianos para una simulación y así He intentado lo siguiente:

import System.Random 

seed = 10101 
gen = mkStdGen seed 

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

Esto es sólo el algoritmo de Box-Muller - r1 dados, uniforme r2 al azar números en el intervalo [0,1] devuelve un número aleatorio gaussiano.

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

así que utiliza esta función cada vez que normals que necesitaba mi lista de números aleatorios.

El problema con eso debe ser evidente: genera siempre la misma secuencia porque siempre estoy usando la misma semilla. No obtengo nuevas secuencias, solo obtengo los primeros n valores de la secuencia todo el tiempo.

Lo que estaba fingiendo claramente fue que, cuando escribo:

x = normal 10 
y = normal 50 

Tendría x ser los primeros 10 valores de map (boxMuller 0 1) $ pairs $ randoms g e y ser los siguientes 50 valores de esta lista, y así en.

Por supuesto, esto es imposible, porque una función siempre debe devolver los mismos valores con la misma entrada. ¿Cómo escapo de esta trampa?

Answer 1

Creo que hacer los cálculos que requieren números aleatorios dentro de una mónada que abstrae el generador es lo más limpio. Aquí está cómo se vería ese código:

Vamos a poner la instancia de StdGen en una mónada de estado, luego proporcionaremos algo de azúcar sobre el método de obtención y configuración de la mónada de estado para darnos números aleatorios.

En primer lugar, cargar los módulos:

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

(Normalmente Lo que probablemente no especificar las importaciones, pero esto hace que sea fácil de entender que cada función está viniendo.)

A continuación, vamos a definir nuestra mónada de "cómputos que requieren números aleatorios"; en este caso, un alias para State StdGen llamado R. (Porque "Aleatorio" y "Rand" ya significan algo más).)

type R a = State StdGen a 

La forma en que funciona R es: definir un cálculo que requiere una serie de números al azar (el monádico "acción"), y luego se "ejecuta" con runRandom:

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

Esto requiere una acción y una semilla, y devuelve los resultados de la acción. Al igual que con el evalState, runReader, etc.

Ahora solo necesitamos azúcar en la mónada estatal. Usamos get para obtener StdGen, y usamos put para instalar un nuevo estado. Esto significa que, para obtener un número aleatorio, escribiríamos:

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

obtenemos el estado actual del generador de números aleatorios, lo utilizan para obtener un nuevo número aleatorio y un nuevo generador, guardar el número al azar, instale el nuevo estado del generador, y devuelve el número aleatorio.

Ésta es una "acción" que se pueden ejecutar con runRandom, por lo que vamos a intentarlo:

ghci> runRandom rand 42 
0.11040701265689151       
it :: Double  

Ésta es una función pura, así que si lo haces otra vez con los mismos argumentos, obtendrá la mismo resultado. La impureza permanece dentro de la "acción" que pasas a runRandom.

De todos modos, su función quiere pares de números aleatorios, por lo que vamos a escribir una acción para producir un par de números aleatorios:

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

Ejecutar este con runRandom, y verá dos números diferentes en la pareja, como era de esperar. Pero note que no tuvo que proporcionar "rand" con una discusión; eso es porque las funciones son puras, pero "rand" es una acción que no necesita ser pura.

Ahora podemos implementar su función normals. boxMuller es como se ha definido lo anterior, me acaba de agregar una firma tipo para que pueda entender lo que está pasando sin tener que leer toda la función:

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

Con todas las funciones de ayuda/acciones llevadas a cabo, por fin podemos escribir normals , una acción de 0 argumentos que devuelve un (generada con pereza-) lista infinita de dobles distribuido normalmente:

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

usted podría también escribir esta concisa menos si quieres:

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

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

repeat() da a la acción monádica una corriente infinita de nada para invocar la función con (y es lo que hace que el resultado de las normales sea infinitamente largo). Inicialmente había escrito [1..], pero lo reescribí para eliminar el 1 sin sentido del texto del programa. No estamos operando en enteros, solo queremos una lista infinita.

De todos modos, eso es todo.Para utilizar esto en un programa real, que acaba de hacer su trabajo que requiere normales al azar dentro de una acción 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 

Las técnicas habituales de acciones monádicos de programación se aplican; mantenga tan poco de su aplicación en R como sea posible, y las cosas no serán demasiado complicadas.

Ah, y en otra cosa: la pereza puede "fugas" fuera del límite de la mónada limpiamente. Esto significa que puede escribir:

take 10 $ runRandom normals 42 

y funcionará.

Answer 2

La lista que obtiene de randoms es infinita, y cuando usa el prefijo finito, no necesita tirar la cola infinita. Puede pasar los números aleatorios como un parámetro adicional y devolver los números aleatorios no utilizados como un resultado adicional, o puede estacionar una secuencia infinita de números aleatorios en una mónada de estado.

Se produce un problema similar para compiladores y otros códigos que desean un suministro de símbolos únicos. Esto es solo un verdadero dolor en Haskell, porque estás enhebrando el estado (del generador de números aleatorios o del generador de símbolos) en todo el código.

He hecho algoritmos aleatorizados tanto con parámetros explícitos como con una mónada, y ninguno de los dos es realmente satisfactorio. Si asimilas mónadas, probablemente tenga una ligera recomendación para usar una mónada de estado que contenga la lista infinita de números aleatorios que aún no se han usado.

Answer 3

También podría eludir el problema utilizando newStdGen y luego obtendrá una semilla diferente (prácticamente) cada vez.