A partir de la nueva versión del estándar de C++ conocida como C++11 (o en su última actualización C++14 http://isocpp.org/blog/2014/08/we-have-cpp14) la librería estándar aporta un buen número de clases y tipos relacionados con la generación de números aleatorios y sus distribuciones. Esta es una novedad de interés para la biología computacional puesto que gran cantidad de algoritmos y simulaciones se apoyan en el modelo de Monte Carlo cuyo pilar fundamental es la generación de números aleatorios. Antes de esta inclusión era el propio desarrollador quién debía encargarse de construir su propia librería de funciones aleatorias o importar alguna de distribución libre. Algunas librerías de este tipo eran las antiguas genlib.c y randlib.c, o, más recientemente, ALGLIB http://www.alglib.net/ y la librería científica de la GNU (GSL) http://www.gnu.org/software/gsl/ que incluían, además de otras muchas funciones matemáticas, también las relacionadas con la generación de números aleatorios y distribuciones de probabilidad.
Sin embargo, la nueva librería de números aleatorios incluida en C++11 además de algunas de las más importantes distribuciones, aporta también algunos de los mejores y más modernos generadores de aleatoriedad. Sin entrar en detalles que nos llevarían más allá de la intención de este post, podemos aclarar que en la generación de variables aleatorias distinguimos dos aspectos fundamentales. Primero, necesitamos un mecanismo generador de aleatoriedad (random engine) es decir una fuente que genere patrones "impredecibles" de bits. En segundo lugar, necesitamos funciones que aprovechando la susodicha fuente de aleatoriedad generen números según alguna distribución de probabilidad deseada. De manera un poco confusa, a la combinación de un generador de aleatoriedad y una distribución de probabilidad se le llama generador de números aleatorios (random number generator, RNG).
Debe notarse que el generador de aleatoriedad es en realidad un generador de pseudoaleatoriedad (PRNG) es decir se genera aleatoriedad desde el punto de vista estadístico pero basándose en un modelo causal y determinista. Una simple función matemática que para un valor de entrada (semilla) genera una secuencia (muy larga) de valores de salida. Estos valores deben ser estadísticamente uniformes e independientes. La impredecibilidad radica en que estas funciones generan resultados muy distintos para pequeñas variaciones en los estados iniciales. Por tanto, si se cambia la semilla, el estado inicial del mecanismo generador cambia y la secuencia subsiguiente de números aleatorios será muy distinta a la anterior. La ventaja del uso de procesos pseudoaleatorios es que el desarrollador puede, para una misma semilla, generar exactamente la misma secuencia de números una y otra vez, lo cual es extremadamente útil a la hora de comprobar y depurar código, corregir errores, etc.
La clave para conseguir emular la deseada impredecibilidad de muchos fenómenos estocásticos será conseguir semillas del modo menos predecible posible. Una práctica muy común es usar los milisegundos (o nanosegundos) del reloj del procesador. Esta práctica sigue siendo perfectamente válida para la mayoría de aplicaciones en biología computacional si bien ha quedado obsoleta para algunos usos como la seguridad informática (por ejemplo el hackeo del sitio web de tecnología Hacker News se pudo hacer gracias a que usaban números pseudolaeatorios incializados con una semilla basada en el reloj del sistema en milisegundos: https://blog.cloudflare.com/why-randomness-matters/).
Por otro lado, para conseguir generadores realmente aleatorios se pueden usar procesos basados en hardware como alteraciones debidas al ruido electrónico o efectos cuánticos, sin embargo, desde el punto de vista de la simulación, se pierde eficiencia, reproducibilidad y portabilidad. Por tanto es el uso de generadores pseudoaleatorios el que está más extendido, siendo el elemento fundamental de la simulación Monte Carlo y por tanto una parte importantísima de la biología computacional.
El lector interesado puede encontrar un buen tutorial sobre la generación de números aleatorios en C++11 en el siguiente enlace: http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2013/n3551.pdf.
En lo que sigue haré un breve resumen sobre los generadores aleatorios en C++11 y algunos ejemplos de su uso en simulación de biología de poblaciones.
Mecanismos generadores de aleatoriedad en la librería estándar de C++11
Según el tipo de algoritmo que se utilice para generar los patrones aleatorios de bits distinguimos entre diferentes clases de generadores. Por ejemplo, los generadores congruentes lineales (http://en.wikipedia.org/wiki/Linear_congruential_generator) y uno de los más usados actualmente, el algoritmo "Mersenne Twister" http://en.wikipedia.org/wiki/Mersenne_Twister.
Recordad que necesitábamos dos componentes para generar variables aleatorias: un mecanismo generador y una distribución. En C++11 vienen predefinidos varios mecanismos generadores de pseudoaleatoriedad, lo mismo que diferentes distribuciones. Un resumen de todo ello se puede consultar en el siguiente enlace http://en.cppreference.com/w/cpp/numeric/random.
Personalmente suelo utilizar uno de los generadores más comunes como es el std::mt19937 o su versión de 64 bits std::mt19937_64 aunque para casos puntuales donde se van a generar pocos valores, como por ejemplo dentro de una función de muestreo sin reemplazamiento, uso el generador por defecto std::default_random_engine.
El siguiente paso es generar la semilla (seed) que inicializa la serie. Una opción es
auto seed = chrono::high_resolution_clock::now().time_since_epoch().count();
que guarda en seed el número de tics de reloj desde el tiempo Unix o epoch (https://es.wikipedia.org/wiki/Tiempo_Unix). Como este número depende de la precisión del reloj del sistema (microsegundos en mi máquina) podemos querer definirlo con una precisión distinta. Si por ejemplo queremos milisegundos haríamos:
auto seed =
chrono::duration_cast < chrono::milliseconds
O nanosegundos:
chrono::duration_cast
Por supuesto, hay mejores (y más complicadas) formas de generar la semilla. Algunas sugerencias pueden consultarse aquí: http://stackoverflow.com/questions/24334012/best-way-to-seed-mt19937-64-for-monte-carlo-simulations.
Distribuciones de probabilidad en C++11
Existen varios tipos genéricos de distribuciones tanto discretas como continuas. Algunos ejemplos son: uniforme (incluye uniforme entera y real), Bernoulli (incluye Bernoulli, binomial, binomial negativa y geométrica), Poisson (incluye Poisson, exponencial y gamma entre otras) y normal (incluye normal, log-normal y chi-cuadrado entre otras).
Hay detalles a los que vale la pena estar atento, como por ejemplo el rango de los valores generado. Así, la uniforme entera genera valores en el intervalo cerrado [a,b] mientras que la uniforme real lo hace en el semicerrado [a,b).
Algunos ejemplos de generación y uso de números aleatorios
Lo primero que hacemos al comienzo de nuestro programa es definir el mecanismo generador deseado, la semilla y las distribuciones que necesitamos:
mt19937_64 gener; // este es el mecanismo generador que usarán las distribuciones
// la semilla:
auto seed = chrono::high_resolution_clock::now().time_since_epoch().count();
gener.seed(seed); // inicializa el generador con la semilla
// uniforme entre 0 (incluido) y 1 (no incluido):
uniform_real_distribution
// la normal de media m y desviación típica s:
double m(0.0), s(1.0);
normal_distribution
// Poisson de media avemut=5.2
double avemut(5.2);
poisson_distribution< int>
Después se utilizarán unas u otras según sea necesario.
Ensayos Bernoulli: supervivencia y fertilidad en una población biológica
Un ensayo Bernoulli con probabilidad p es un experimento con sólo dos resultados posibles donde la probabilidad de uno de ellos es p y la del otro 1 – p. Podemos imaginar muchos procesos biológicos como ensayos Bernoulli. Por ejemplo, la probabilidad de sobrevivir dado un valor de eficacia w o la probabilidad de ser progenitor dado un valor de fertilidad f etc etc. Consideremos la eficacia de un individuo en su aspecto de viabilidad es decir, la probabilidad de sobrevivir hasta alcanzar la madurez sexual. Entonces, la supervivencia de un individuo X cualquiera con valor de eficacia w se puede modelar como una Bernoulli. La probabilidad de que ese individuo sobreviva vendrá dada por
if(X.w > uniform(gener)) ... ; // sobrevive
Nótese que si X.w = 1 el individuo siempre sobrevive (recuerda que el intervalo de la uniforme real era semicerrado) pero si X.w es 0 nunca lo hará.
Inicializar los individuos de una población con un fenotipo normalmente distribuido
Esta es una asunción muy común en biología. Una población P en la que estamos considerando un carácter cualquiera: color, tamaño, altura, etc... se distribuye de manera normal para ese carácter. Si queremos generar una población inicial de media m y desviación s hacemos:
for(auto indiv=0; indiv< P.size(); ++indiv)
P[indiv].caracter= normal(gener);
Generar los individuos que van a mutar en una población
Dado que la ocurrencia de una mutación es un suceso de muy baja probabilidad se suele modelar mediante la distribución de Poisson. Si la probabilidad de mutación en un genoma haploide es 0.01 y tenemos N individuos diploides entonces esperamos 0.02N nuevas mutaciones en cada generación. Podemos forzar la ocurrencia de ese número exacto de nuevos mutantes en cada generación o, si queremos ser más realistas, dejar que varíe según una Poisson definida con el parámetro 0.02N. De este último modo el número de mutantes oscilará en cada generación aunque el valor medio a lo largo de las generaciones tenderá a 0.02N.
mutantes=poisson(gener);
Estos son sólo algunos ejemplos bastante sencillos de entre la multitud de posibilidades que los generadores aleatorios nos ofrecen para simular procesos estocásticos en biología.
Un último apunte respecto a la semilla es que, si estamos realizando varias réplicas de un proceso biológico relativamente largo, puede ser conveniente iniciar cada réplica con una nueva semilla de modo que actualizamos el generador:
// actualiza el generador si se quieren realizar nuevas simulaciones independientes
gen.seed(chrono::high_resolution_clock::now().time_since_epoch().count());
Contrastar la corrección de un generador de números pseudoaleatorios
Es siempre conveniente comprobar que hemos implementado correctamente nuestros generadores. Para ello podemos comparar valores esperados de cada distribución, como su media y varianza, con los obtenidos en una secuencia dada de números pseudoaleatorios.
La media de una uniforme U(0,1) es 1/2 y la varianza 1/12. Si generamos una cantidad n de valores bajo esta distribución, la media y varianza de esos valores deberían ir aproximándose a los esperados según aumenta n.
Lo mismo podemos hacer con otras distribuciones. Por ejemplo si un individuo tiene viabilidad 0.7 esa es su probabilidad de supervivencia. Si repetimos el proceso n veces esperamos que el número de veces que sobreviva se aproxime a 0.7n y la varianza de los valores de sobrevivir (1) o morir (0) se aproximará a 0.21n. Igualmente, si inicializamos una población de fenotipos utilizando una normal de media m y desviación s, el fenotipo medio y la desviación en la población obtenida deberían aproximarse a m y s según aumenta el tamaño de la población.
Para una información más detallada sobre estos y otros tipos de test, puede consultarse http://www.euclides.dia.uned.es/aurquia/Files/Simulacion%20Teoria%202008_09.pdf
Bibliografia
- The C++ Standard Library. A Tutorial and Reference. Autor: Nicolai Josuttis. Ed. Addison-Wesley 2012.
- Random Number Generation in C++11.Document #: WG21 N3551. url: http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2013/n3551.pdf. Date: 2013-03-12.Revises: None.Project: JTC1.22.32 Programming Language C++.Walter E. Brown.
- Matsumoto, M.; Nishimura, T. (1998). "Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator". ACM Transactions on Modeling and Computer Simulation 8 (1): 3–30.
- Urquia: "Simulacion, Texto base de teoria", 200.url: http://www.euclides.dia.uned.es/aurquia/Files/Simulacion%20Teoria%202008_09.pdf
