triangulo - Muestreo de secuencias de números aleatorios en Haskell

numeros complejos en haskell

Question

Necesito pequeñas listas de números aleatorios gaussianos para una simulación y por eso probé 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) 

Este es solo el algoritmo de Box-Muller: dado los números aleatorios uniformes r1, r2 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 utilicé esta función normals cada vez que necesitaba mi lista de números aleatorios.

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

Lo que pretendía claramente era que, cuando escribía:

x = normal 10 
y = normal 50

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

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

Answer 1

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

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

Primero, cargue los módulos:

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

(Normalmente, es probable que no especifique las importaciones, pero esto hace que sea fácil entender de dónde proviene cada función).

Luego definiremos nuestra mónada "cálculos que requieren números aleatorios"; en este caso, un alias para State StdGen llamado R (Porque "Aleatorio" y "Rand" ya significan otra cosa.)

type R a = State StdGen a

La forma en que R funciona es: usted define un cálculo que requiere un flujo de números aleatorios (la "acción" monádica), y luego lo "ejecuta" con runRandom :

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

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

Ahora solo necesitamos azúcar alrededor de la mónada estatal. Usamos get para obtener el 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 usamos para obtener un nuevo número aleatorio y un nuevo generador, guardamos el número aleatorio, instalamos el nuevo estado del generador y devolvemos el número aleatorio.

Esta es una "acción" que se puede ejecutar con runRandom, así que intentémoslo:

ghci> runRandom rand 42
0.11040701265689151                           
it :: Double     

Esta es una función pura, así que si la ejecutas de nuevo con los mismos argumentos, obtendrás el 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, así que escribamos una acción para producir un par de números aleatorios:

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

Ejecuta esto con runRandom, y verás dos números diferentes en el par, como es de esperar. Pero tenga en cuenta que no tuvo que proporcionar a "rand" un argumento; 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 lo boxMuller anteriormente, acabo de agregar una firma de 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 / acciones de ayuda implementadas, finalmente podemos escribir normals , una acción de 0 argumentos que devuelve una lista infinita (generada perezosamente) de dobles distribuidos normalmente:

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

También puedes escribir esto de manera menos concisa si quieres:

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

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

repeat () le da a la acción monádica un flujo infinito 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 significado del texto del programa. No estamos operando en enteros, solo queremos una lista infinita.

De todos modos, eso es todo. Para usar esto en un programa real, solo hace su trabajo requiriendo normales aleatorias 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

Se aplican las técnicas habituales para programar acciones monádicas; mantenga la menor cantidad posible de su aplicación en R , y las cosas no serán demasiado complicadas.

Ah, y por otro lado: la pereza puede "filtrarse" fuera del límite de la mónada limpiamente. Esto significa que puedes escribir:

take 10 $ runRandom normals 42

y funcionará.

Answer 2

La lista que obtiene de los 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 estatal.

Un problema similar ocurre con los compiladores y otros códigos que desean un suministro de símbolos únicos. Esto es solo una molestia real en Haskell, porque está en estado de hebra (del generador de números aleatorios o del generador de símbolos) en todo el código.

He hecho algoritmos aleatorios con parámetros explícitos y con una mónada, y ninguno de los dos es realmente satisfactorio. Si no te gustan las mónadas, probablemente te recomiendo que uses una mónada estatal que contenga la lista infinita de números aleatorios que aún no se han utilizado.

Answer 3

También puede evitar el problema usando newStdGen y luego obtendrá una semilla diferente (virtualmente) cada vez.