Distribuciones discretas con mónadas

Table of Contents

El modelo

Como parte de una serie de ejemplos sobre uso de mónadas, he escrito un poco de código para modelar distribuciones discretas usando mónadas. Por un lado, usa un generador congruencial lineal para generar números aleatorios; y por otro, usa la mónada State para pasar una semilla aleatoria de una función a otra que me permita seguir generado números aleatorios. Por último, aporta un método que deriva Show para probar las distribuciones y dibujar un histograma de cualquiera de ella.

Componiendo distribuciones

Lo más útil de esta idea es el poder generar unas distribuciones a partir de otras. La primera que intentamos es una uniforme discreta (un dado de n caras) usando una semilla inicial. En el siguiente código se implementa el generador congruencial.

  dice :: Int -> Distribution Int
  dice n = state (\s -> (s `mod` n + 1, 16807*s `mod` 2147483647))

Vemos que funciona como una distribución uniforme.

  >>> dice 6
  
  1:     ################
  2:     ################
  3:     ################
  4:     ################
  5:     ################
  6:     ################

Y desde ella generar fácilmente otras usando funciones que compongan distribuciones. Un ejemplo es usar (⊕) = liftM2 (+) para sumar dados.

>>> dice 6 ⊕ dice 6

2:   #####
3:   ##########
4:   ###############
5:   ####################
6:   ##########################
7:   ##############################
8:   #########################
9:   ####################
10:  ###############
11:  ##########
12:  #####

Otras distribuciones

Si seguimos componiendo usando la estructura de mónada, podemos crear otras distribuciones simples como la distribución de Bernoulli y la distribución binomial.

  bernoulli :: Double -> Distribution Int
  bernoulli p = do
    sample <- dice 1000000
    if (fromIntegral sample / 1000000.0 < p)
      then return 1
      else return 0

  binomial :: Int -> Double -> Distribution Int
  binomial k p = sum <$> replicateM k (bernoulli p)

Lo interesante de este código es que dejamos a la estructura de mónada encargarse internamente de el paso de la semilla de aleatoriedad y la construcción de distribuciones complejas puede hacerse composicionalmente.

El código

El siguiente código es una primera implementación de este post en Haskell.

  {-# LANGUAGE FlexibleInstances #-}
  {-# LANGUAGE TypeSynonymInstances #-}
  {-
  En este archivo vamos a usar mónadas para definir distribuciones
  discretas de probabilidad y aplicar operaciones algebraicas sobre ellas.
  -}
  import Control.Monad.State
  
  -- Generación aleatoria
  -- Para generar números pseudoaleatorios usaremos LCGs. La idea es tener
  -- un dado que nos dé una distribución de probabilidad uniforme dada una
  -- semilla y nos devuelva el resultado de la tirada y una nueva semilla
  -- aleatoria. Buscamos que un dado de seis caras sea, por ejemplo:
  --
  --   dice 6 :: Seed -> (Int, Seed)
  --
  -- Si quisiéramos tirar dos dados, tendríamos que tomar la semilla resultante
  -- del primer lanzamiento y pasarla al segundo; algo así:
  --
  --   let (a,newseed) = dice 6 seed
  --   let (b,_)       = dice 6 newseed
  --   print [a,b]
  --
  -- Pero esto se hace demasiado complejo. La semilla, en el fondo, es un
  -- estado, así que podemos modelarla con la mónada State. Cada lanzamiento
  -- será de la forma:
  --
  --   State Seed a   ===   Seed -> (a, Seed)
  --
  -- Luego podemos llamar a la distribución: Distribution a = State Seed a, y
  -- trabajar con ella usando las funciones normales de mónadas.
  type Seed = Int
  type Distribution = State Seed
  
  
  
  -- Nuestra primera distribución es un dado de "n" lados que usa internamente un
  -- generador de números aleatorios.
  dice :: Int -> Distribution Int
  dice n = state (\s -> (s `mod` n + 1, 16807*s `mod` 2147483647))
  
  -- Una moneda es un dado de dos caras
  coin :: Distribution Int
  coin = dice 2
    
  -- Estas funciones pueden ser llamadas con la mónada estado, dada una
  -- semilla inicial, devuelven el resultado y la nueva semilla:
  --
  -- λ> runState (dice 6) 1
  -- (2,16807)
  -- λ> runState (dice 6) 16807
  -- (2,282475249)
  --
  -- El usar composición con mónadas nos ahorraba controlar los errores
  -- en el primer caso, aquí nos ahorra controlar el cambio de semilla,
  -- por ejemplo: para lanzar dos dados y hacer que la semilla se pase
  -- internamente.
  twodices' :: Distribution Int
  twodices' = do
    a <- dice 6
    b <- dice 6
    return (a+b)

  () :: Distribution Int -> Distribution Int -> Distribution Int
  () = liftM2 (+)
  () :: Distribution Int -> Distribution Int -> Distribution Int
  () = liftM2 (*)
  
  twodices :: Distribution Int
  twodices = dice 6  dice 6
  
  -- Igual que hago esto, podría hacer:
  --
  --   foldr (⊕) (return 0) [dice 6,dice 6,dice 6]
  --   foldr (⊕) (return 0) (replicate 10 (dice 6))
  --
  -- Que da un resultado que se aproxima a una distribución normal.
  
  -- Ahora, desde ella, podemos crear otras distribuciones. La distribución de
  -- bernoulli sería la de una moneda trucada donde una cara tiene probabilidad
  -- p y la otra tiene probabilidad (1-p).
  bernoulli :: Double -> Distribution Int
  bernoulli p = do
    sample <- dice 1000000
    if (fromIntegral sample / 1000000.0 < p)
      then return 1
      else return 0

  -- La distribución binomial es la suma de k distribuciones de Bernoulli
  binomial :: Int -> Double -> Distribution Int
  binomial k p = sum <$> replicateM k (bernoulli p)
  
  -- La distribución constante y otra forma de escribir la distribución
  -- binomial, de manera algebraica.
  constant :: Int -> Distribution Int
  constant n = return n
  
  binomial' :: Int -> Double -> Distribution Int
  binomial' k p = foldr () (constant 0) (replicate k (bernoulli p))
  
  
  
  -- Muestra la distribución. Los detalles de implementación no son interesantes.
  -- Hemos usado  TypeSynonymInstances para simplificar el proceso de sobrecargar
  -- la instancia de Show y poder dibujar directamente por la pantalla las
  -- demostraciones.
  instance Show (State Seed Int) where
    show = showdist
  
  showdist :: Distribution Int -> String
  showdist dist = unlines $ map counter [minimum samples..maximum samples]
    where samples = fst $ runState (replicateM 50000 dist) 1
          counter n = show n ++ ":\t " ++ replicate ((count n samples) `div` (3000 `div` range)) '#'
          range = maximum samples - minimum samples + 1
  
  count :: Eq a => a -> [a] -> Int
  count x = length . filter (x==)
  
  
  main :: IO ()
  main = return ()

Last edited 2020-04-25 12:12:11 by Mario Román (code).