Números pseudoaleatorios: C++ Linux

Pseudoaleatory numbers: C++ Linux

Cuando tenemos la necesidad de simular algún fenómeno estocástico (basado en probabilidades) posiblemente habrá que recurrir a los números pseudoaleatorios. Estos números son generados por rutinas que intentan reproducir la propiedad estadística de que cualquiera de ellos tiene misma probabilidad de ocurrencia. Lo de pseudoaleatorio proviene del hecho de que la serie siempre será igual cuando se utiliza el mismo número generador ó semilla (seed en inglés). No obstante, uno de los principales obstáculos que se les presenta a los programadores es que las rutinas deben producir una serie lo suficientemente larga para los propósitos que se requieran, es decir, que no sean cíclicos. Estos números pseudoaleatorios, al igual que las probabilidades, están comprendidos entre 0 y 1. Sin embargo, para algunos propósitos especiales podemos requerir que ellos presenten, adicionalmente, la propiedad de tener media 0 (cero) y varianza 1. Supongamos que queremos generar una secuencia de lluvia. Primero, debemos escoger un modelo estadístico para ello. Estos módelos estadísticos se basan en cadenas de Markov donde la probabilidad de ocurrencia de la precipitación depende de si los días anteriores fueron lluviosos ó no. Sin entrar en las dificultades propias de como encontrarlas diremos, por ahora, que las probabilidades de día húmedo (lluvioso) corresponden a un número entre 0 y 1 y sus valores serán más elevados mientras más lluvioso sea el mes del año al cual hagamos referencia. Por tanto, los meses secos tendrán probabilidades de día húmedo muy bajas y el sentido común nos señala que la mayoría de los números pseudoaleatorios generados en estos casos serán mayores. En consecuencia, las rutinas para generar las secuencias deben comparar el número pseudoaleatorio producido con la probabilidad mensual de día húmedo y si es mayor que ésta entonces el día será seco (se le asigna un cero). En caso contrario el día será húmedo (se le asigna un 1).

Por otra parte, una vez que tenemos que un día es húmedo (unos en la secuencia), se requiere de un modelo estadístico-matemático que nos permita asignarle la magnitud de la cantidad llovida para ese día. Las distribuciones mensuales de lluvia siguen patrones asimétricos (no son distribuciones normales ó gaussianas) que son caracterizados por parámetros estadísticos y se requieren números pseudoaleatorios que presentan las propiedades anteriormente referidas: media cero y varianza uno.

A continuación, incluyo un pequeño programa escrito en C++ cuya rutina para generar los números pseudoaleatorios anteriormente mencionados, la he utilizado en programas para obtener series sintéticas de precipitación.

#include <iostream>
#include <iomanip>
#include <fstream>

using namespace std;

#include <math.h>

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1 + (IM - 1) / NTAB)
#define EPS 1.2e-7
#define RNMX (1.0 - EPS)

double ran1(long *idum);

double gasdev(long *idum);

int main() {

 ofstream label1 ("salida.out");

 long semilla; 

 long *seed = &semilla;

 semilla = -2000000;

 int i;

 cout.setf(ios::fixed);
 cout.setf(ios::right);
 cout.precision(7);

 label1.setf(ios::fixed);
 label1.setf(ios::right);
 label1.precision(7);

 for (i=0; i < 65000; i++){

  cout << setw(10) << ran1(seed) << setw(12) << gasdev(seed) << endl;
  label1 << setw(10) << ran1(seed) << setw(12) << gasdev(seed) << endl;

 } 

 cout << endl;

 return 0;

}

double ran1(long *idum) {

 int j;

 long k;

 static long iy = 0;

 static long iv[NTAB];

 double temp;

 if (*idum <= 0 || !iy) {

  if (-(*idum) < 1) *idum = 1;

   else *idum = -(*idum);

   for (j = NTAB + 7; j>=0; j--) {

    k = (*idum) / IQ;

    *idum = IA * (*idum - k*IQ) - IR*k;

    if (*idum < 0) *idum += IM;

     if (j< NTAB) iv[j] = *idum;

   }

   iy = iv[0];

 }

 k = (*idum) / IQ;

 *idum = IA * (*idum - k * IQ) - IR * k;

 if (*idum < 0) *idum += IM;

 j = iy / NDIV;

 iy = iv[j];

 iv[j] = *idum;

 if ((temp = AM*iy) > RNMX) return RNMX;

  else return temp;

}

double gasdev(long *idum) {

 double ran1(long *idum);

 static int iset =0;

 static double gset;

 double fac, rsq, v1, v2;

 if (*idum < 0) iset = 0;

 if (iset == 0) {

  do {

   v1 = 2.0 * ran1(idum) - 1.0;

   v2 = 2.0 * ran1(idum) - 1.0;

   rsq = v1*v1 + v2*v2;

  } while (rsq >= 1.0 || rsq == 0);

  fac = sqrt (-2.0*log(rsq)/rsq);

  gset = v1 * fac;

  iset = 1;

 return v2*fac;

 } else {

  iset = 0;

  return gset;

 }

}

Después de desempaquetado el programa, nos movemos hacia la carpeta ALEATORIO y allí compilamos con:

g++ aleatorio.c++ -o aleatorio

La ejecución del mismo con ./aleatorio produce una serie pareada (65.000 x 2 = 130.000) de números pseudoaleatorios en salida.out (semilla -2.000.000) cuyas propiedades estadísticas pueden ser verificadas en una hoja de cálculo de Openoffice ó en Excel. La primera serie (primera columna) tiene media 0,5 y la segunda serie (segunda columna) media 0 y varianza 1; tal como se deseaba. Aunque en este caso no lo tomé en cuenta, siempre es conveniente que las semillas se almacenen en archivos externos a fin de que no tengamos que modificar el programa para obtener una nueva serie con otra semilla.

About these ads
Esta entrada fue publicada en Código C++. Guarda el enlace permanente.

Deja un comentario

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s