sequenze di campionamento di numeri casuali a Haskell

Question

ho bisogno di piccole liste di numeri casuali gaussiana per una simulazione e così ho provato quanto segue:

import System.Random 

seed = 10101 
gen = mkStdGen seed 

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

Questo è solo l'algoritmo di Box-Muller - dato r1, r2 uniforme casuale numeri nell'intervallo [0,1] restituisce un numero casuale 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) 

Così ho usato questo normals funzione ogni volta che avevo bisogno di mia lista di numeri casuali.

Il problema deve essere evidente: genera sempre la stessa sequenza perché sto usando sempre lo stesso seme! Non sto ottenendo nuove sequenze, sto solo ottenendo i primi n valori della sequenza tutto il tempo.

Quello che mi fingevo chiaramente era che, quando digitate:

x = normal 10 
y = normal 50 

avrei x essere i primi 10 valori di map (boxMuller 0 1) $ pairs $ randoms g e y essere i 50 valori prossimi in questo elenco, e così sopra.

Ovviamente questo è impossibile, perché una funzione deve sempre restituire gli stessi valori in base allo stesso input. Come posso sfuggire a questa trappola?

Answer 1

Penso che fare i calcoli che richiedono numeri casuali all'interno di una monade che astrae il generatore è la cosa più pulita. Ecco come dovrebbe apparire il codice:

Metteremo l'istanza StdGen in una monade di stato, quindi forniremo dello zucchero sul metodo get e set del monad di stato per fornirci numeri casuali.

In primo luogo, caricare i moduli:

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

(Normalmente avrei probabilmente non specificare le importazioni, ma questo lo rende facile da capire dove ogni funzione è venuta da.)

Poi ci definirà i nostri calcoli "che richiedono numeri casuali"; in questo caso, un alias per State StdGen chiamato R. (Perché "Random" e "Rand" significano già qualcos'altro.)

type R a = State StdGen a 

Il modo in cui R funziona è: si definisce un calcolo che richiede un flusso di numeri casuali (l ' "azione" monadica), e poi si "esegue" con runRandom:

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

Ciò richiede un'azione e un seme e restituisce i risultati dell'azione. Proprio come i soliti evalState, runReader, ecc.

Ora abbiamo solo bisogno di zucchero attorno alla monade di Stato. Usiamo get per ottenere StdGen e usiamo put per installare un nuovo stato. Ciò significa, per ottenere un numero casuale, scriveremo:

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

otteniamo lo stato corrente del generatore di numeri casuali, utilizzarlo per ottenere un nuovo numero casuale e un nuovo generatore, salvare il numero casuale, installare il nuovo stato del generatore e restituisce il numero casuale.

Si tratta di un "" azione" che può essere eseguito con runRandom, quindi proviamo:

ghci> runRandom rand 42 
0.11040701265689151       
it :: Double  

Questa è una funzione di pura, quindi se lo si esegue di nuovo con gli stessi argomenti, si otterrà il stesso risultato. L'impurità rimane all'interno dell'azione che si passa a runRandom.

In ogni caso, la funzione vuole coppie di numeri casuali, quindi cerchiamo di scrivere un'azione per produrre un due di numeri casuali:

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

Esegui questo con runRandom, e vedrete due numeri diversi in la coppia, come ti aspetteresti. Ma notate che non avete dovuto fornire "rand" con un argomento; questo perché le funzioni sono pure, ma "rand" è un'azione, che non deve essere pura.

Ora possiamo implementare la tua funzione normals. boxMuller è come definito sopra, ho solo aggiunto una firma di tipo così posso capire quello che sta succedendo, senza dover leggere l'intera funzione:

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

Con tutte le funzioni di supporto/azioni attuate, possiamo finalmente scrivere normals , un'azione di 0 argomenti che restituisce un (pigramente-generata) lista infinita di normalmente distribuito Doppio:

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

si potrebbe anche scrivere questo meno conciso se si desidera:

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

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

repeat() conferisce all'azione monadica un flusso infinito di nulla per invocare la funzione con (ed è ciò che rende infinitamente lungo il risultato delle normali). Avevo inizialmente scritto [1..], ma l'ho riscritto per rimuovere l'inutile 1 dal testo del programma. Non stiamo operando su interi, vogliamo solo una lista infinita.

In ogni caso, è tutto.Per utilizzare questo in un vero e proprio programma, basta fare il vostro lavoro richiede normali casuali all'interno di un'azione 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 

Le tecniche usuali per la programmazione delle azioni monadici applicano; mantieni il più possibile la tua applicazione nel R e le cose non saranno troppo confuse.

Oh, e su altra cosa: la pigrizia può "trapelare" al di fuori del confine della monade in modo pulito. Ciò significa che è possibile scrivere:

take 10 $ runRandom normals 42 

e funzionerà.

Answer 2

L'elenco che si ottiene da randoms è infinito e quando si utilizza il prefisso finito, non è necessario eliminare la coda infinita. È possibile passare i numeri casuali come parametro aggiuntivo e restituire i numeri casuali non utilizzati come risultato aggiuntivo, oppure è possibile parcheggiare una sequenza infinita di numeri casuali in una monade di stato.

Un problema simile si verifica per i compilatori e altri codici che desiderano fornire simboli univoci. Questo è solo un vero dolore per Haskell, perché si sta eseguendo il threading dello stato (del generatore di numeri casuali o del generatore di simboli) in tutto il codice.

Ho eseguito algoritmi randomizzati sia con parametri espliciti che con una monade, e nessuno dei due è davvero soddisfacente. Se batti monadi, probabilmente ho una leggera raccomandazione di usare una monade di stato contenente l'infinita lista di numeri casuali che non sono ancora stati usati.

Answer 3

Si potrebbe anche aggirare il problema utilizzando newStdGen e quindi si otterrà un seme diverso (praticamente) ogni volta.