Crear y controlar una secuencia de números aleatorios

La traducción de esta página está obsoleta. Haga clic aquí para ver la última versión en inglés.

La clase RandStream le permite crear una secuencia de números aleatorios. Esto resulta útil por varias razones:

Puede generar valores aleatorios sin que afecte al estado de la secuencia global.
Puede separar las fuentes de aleatoriedad en una simulación.
Puede utilizar un generador que se haya configurando de manera diferente que al que utiliza el software de MATLAB^® al iniciarse.

Con el objeto RandStream, puede crear su propia secuencia, establecer los valores de las propiedades y utilizarla para generar números aleatorios. Puede controlar la secuencia que cree de la misma manera que controla la secuencia global. Incluso puede sustituir la secuencia global por la secuencia que cree.

Para crear una secuencia, utilice la función RandStream.

myStream = RandStream('mlfg6331_64');
rand(myStream,1,5)

ans =
    0.6986    0.7413    0.4239    0.6914    0.7255

La secuencia aleatoria myStream actúa de forma independiente a la secuencia global. Si llama a las funciones rand, randn, randi y randperm con myStream como el primer argumento, extraen de las secuencias que ha creado. Si llama a rand, randn, randi y randperm sin myStream, extraen de la secuencia global.

Puede hacer que myStream sea la secuencia global con el método RandStream.setGlobalStream.

RandStream.setGlobalStream(myStream)
RandStream.getGlobalStream

ans = 

mlfg6331_64 random stream (current global stream)
             Seed: 0
  NormalTransform: Ziggurat

RandStream.getGlobalStream == myStream

ans =
     1

Subsecuencias

Puede utilizar subsecuencias para obtener resultados diferentes que sean independientes estadísticamente de una secuencia. A diferencia de las semillas, donde las ubicaciones en la secuencia de números aleatorios no se conocen exactamente, el espacio entre las subsecuencias es conocido, por lo que cualquier posibilidad de superposición se puede eliminar. En resumen, las subsecuencias son una forma más controlada de hacer muchas de las mismas cosas para las que tradicionalmente se han usado las semillas. Las subsecuencias también son una solución más ligera que las secuencias paralelas.

Las subsecuencias proporcionan una manera más rápida y sencilla de garantizar que obtiene resultados diferentes para el mismo código en momentos diferentes. Para utilizar la propiedad Substream de un objeto RandStream, cree una secuencia con un generador que admita subsecuencias. Consulte Elegir un generador de números aleatorios para obtener una lista de algoritmos de generación compatibles con las subsecuencias y sus propiedades. Por ejemplo, genere varios números aleatorios en un bucle.

myStream = RandStream('mlfg6331_64');
RandStream.setGlobalStream(myStream)
for i = 1:5
    myStream.Substream = i;
    z = rand(1,i)
end

z =
    0.6986

z =
    0.9230    0.2489

z =
    0.0261    0.2530    0.0737

z =
    0.3220    0.7405    0.1983    0.1052

z =
    0.2067    0.2417    0.9777    0.5970    0.4187

En otro bucle, puede generar valores aleatorios que sean independientes del primer conjunto de 5 iteraciones.

for i = 6:10
    myStream.Substream = i;
    z = rand(1,11-i)
end

z =
    0.2650    0.8229    0.2479    0.0247    0.4581

z =
    0.3963    0.7445    0.7734    0.9113

z =
    0.2758    0.3662    0.7979

z =
    0.6814    0.5150

z =
    0.5247

Las subsecuencias son útiles en cálculos de serie. Las subsecuencias pueden recrear una simulación completa o parte de ella regresando a un punto de control particular de una secuencia. Por ejemplo, puede regresar a la 6.ª subsecuencia en el bucle. El resultado contiene los mismos valores que el 6.º resultado anterior.

myStream.Substream = 6;
z = rand(1,5)

z =
    0.2650    0.8229    0.2479    0.0247    0.4581

Elegir un generador de números aleatorios

MATLAB ofrece varias opciones de algoritmos de generación. La tabla resume las propiedades esenciales de los algoritmos de generación y las palabras clave que se usan para crearlos. Para recuperar una lista de todos los algoritmos de generación disponibles, utilice el método RandStream.list.

Palabra clave	Generador	Soporte de secuencias y subsecuencias múltiples	Periodo aproximado con plena precisión
`mt19937ar`	Mersenne twister (se utiliza de forma predeterminada al iniciar MATLAB)	No	2¹⁹⁹³⁷-1
`dsfmt19937`	Mersenne twister rápido orientado a SIMD	No	2¹⁹⁹³⁷-1
`mcg16807`	Generador congruencial multiplicativo	No	2³¹-2
`mlfg6331_64`	Generador de Fibonacci retardado multiplicativo	Sí	2¹²⁴ (2⁵¹ secuencias de longitud 2⁷²)
`mrg32k3a`	Generador combinado múltiple recursivo	Sí	2¹⁹¹ (2⁶³ secuencias de longitud 2¹²⁷)
`philox4x32_10`	Generador Philox 4x32 con 10 rondas	Sí	2¹⁹³ (2⁶⁴ secuencias de longitud 2¹²⁹)
`threefry4x64_20`	Generador Threefry 4x64 con 20 rondas	Sí	2⁵¹⁴ (2²⁵⁶ secuencias de longitud 2²⁵⁸)
`shr3cong`	Generador de registro de desplazamientos sumado a un generador lineal congruencial	No	2⁶⁴
`swb2712`	Generador modificado de resta con préstamo	No	2¹⁴⁹²

Los generadores mcg16807, shr3cong y swb2712 son compatibles con versiones anteriores de MATLAB. mt19937ar y dsfmt19937 están diseñados principalmente para aplicaciones secuenciales. Los generadores restantes proporcionan soporte específico para la generación de números aleatorios paralela.

Según la aplicación, puede que ciertos generadores sean más rápidos o devuelvan valores con mayor precisión. Todos los generadores de números seudoaleatorios se basan en algoritmos deterministas y todos los generadores pasarán una prueba suficientemente específica de aleatoriedad. Una forma de comprobar los resultados de la simulación Montecarlo es volver a ejecutar la simulación con dos o más algoritmos de generación y la oferta de generadores que le ofrece MATLAB le permite hacerlo. Aunque al usar generadores diferentes es poco probable que los resultados difieran en más de un error de muestra de Montecarlo, en la literatura hay ejemplos en los que este tipo de validación ha dado lugar a fallos en un algoritmo de generación concreto. (Para ver un ejemplo, consulte [13]).

Algoritmos de generación

mt19937ar: El Mersenne twister, como se describe en [11], tiene periodo $2^{19937} - 1$ y cada valor U(0,1) se crea con números enteros de 32 bits. Los posibles valores son múltiplos de $2^{- 53}$ en el intervalo (0, 1). Este generador no es compatible con secuencias ni subsecuencias múltiples. El algoritmo randn que se utiliza de forma predeterminada para las secuencias mt19937ar es el algoritmo zigurat [7], pero con el generador mt19937ar debajo.
Nota
Este generador es idéntico al que utiliza la función rand a partir de la Versión 7 de MATLAB, que se activa con rand('twister',s).
dsfmt19937: El Mersenne twister rápido orientado a SIMD de precisión doble, como se describe en [12], es una implementación más rápida del algoritmo Mersenne twister. El periodo es $2^{19937} - 1$ y los posibles valores son múltiplos de $2^{- 52}$ en el intervalo (0, 1). De forma nativa, el generador produce valores de precisión doble en [1, 2), que se transforman para crear valores U(0, 1). Este generador no es compatible con secuencias ni subsecuencias múltiples.
mcg16807: Un generador congruencial multiplicativo de 32 bits, como se describe en [14], con multiplicador $a = 7^{5}$ , módulo $m = 2^{31} - 1$ . Este generador tiene un periodo de $2^{31} - 2$ y no admite secuencias ni subsecuencias múltiples. Cada valor U(0, 1) se crea con un número entero de 32 bits a partir del generador; los valores posibles son todos múltiplos de ${(2^{31} - 1)}^{- 1}$ estrictamente dentro del intervalo (0, 1). Para las secuencias de mcg16807, el algoritmo que utiliza randn es el algoritmo polar (descrito en [1]).
Nota
Este generador es idéntico al que usan las funciones rand y randn a partir de la Versión 4 de MATLAB, que se activa con rand('seed',s) o con randn('seed',s).
mlfg6331_64: Un generador de Fibonacci retardado multiplicativo de 64 bits, como se describe en [10], con retardos $l = 63$ , $k = 31$ . Este generador es similar al MLFG implementado en el paquete SPRNG. Tiene un periodo de aproximadamente $2^{124}$ . Admite hasta $2^{61}$ secuencias paralelas, a través de la parametrización, y $2^{51}$ subsecuencias, cada una de longitud $2^{72}$ . Cada valor U(0,1) se crea con el número entero de 64 bits del generador; los valores posibles son todos múltiplos de $2^{- 64}$ estrictamente dentro del intervalo (0, 1). El algoritmo randn que se utiliza de forma predeterminada para las secuencias mlfg6331_64 es el algoritmo zigurat [7], pero con el generador mlfg6331_64 debajo.
mrg32k3a: Un generador combinado múltiple recursivo de 32 bits, como se describe en [2]. El generador es similar al CMRG implementado en el paquete RngStreams en C. Tiene un periodo de $2^{191}$ y admite hasta $2^{63}$ secuencias paralelas a través de la parametrización, cada una de longitud $2^{127}$ . También admite $2^{51}$ subsecuencias, cada una de una longitud $2^{76}$ . Cada valor U(0, 1) se crea con dos números enteros de 32 bits a partir del generador; los valores posibles son múltiplos de $2^{- 53}$ estrictamente dentro del intervalo (0, 1). El algoritmo randn que se utiliza de forma predeterminada para las secuencias mrg32k3a es el algoritmo zigurat [7], pero con el generador mrg32k3a debajo.
philox4x32_10: Un generador 4 por 32 con 10 rondas, como se describe en [15]. Este generador utiliza una red Feistel y multiplicación de números enteros. El generador está específicamente diseñado para alto rendimiento en sistemas altamente paralelos como las GPU. Tiene un periodo de 2¹⁹³ (2⁶⁴ secuencias de longitud 2¹²⁹).
threefry4x64_20: Un generador 4 por 64 con 20 rondas, como se describe en [15]. Este generador es una adaptación no criptográfica del cifrado de bloque Threefish a partir de la función hash Skein. Tiene un periodo de 2⁵¹⁴ (2²⁵⁶ secuencias de longitud 2²⁵⁸).
shr3cong: El generador SHR3 de Marsaglia de registro de desplazamientos sumado a un generador lineal congruencial con multiplicador $a = 69069$ , sumando $b = 1234567$ y módulo $2^{- 32}$ . SHR3 es un generador de registro de 3 desplazamientos definido como $u = u (I + L^{13}) (I + R^{17}) (I + L^{5})$ , donde $I$ es el operador de identidad, $L$ es el operador de desplazamiento hacia la izquierda y R es el operador de desplazamiento hacia la derecha. El generador combinado (la parte SHR3 se describe en [7]) tiene un periodo de aproximadamente $2^{64}$ . Este generador no es compatible con secuencias ni subsecuencias múltiples. Cada valor U(0,1) se crea con el número entero de 32 bits del generador; los valores posibles son todos múltiplos de $2^{- 32}$ estrictamente dentro del intervalo (0, 1). El algoritmo randn que se utiliza por defecto para las secuencias shr3cong es la forma anterior del algoritmo zigurat [9], pero con el generador shr3cong debajo. Este generador es idéntico al que utiliza la función randn a partir de la Versión 5 de MATLAB, que se activa con randn('state',s).
Nota
El generador SHR3 que se utiliza en [6] (1999) es diferente al que se usa en [7] (2000). MATLAB utiliza la versión más reciente del generador, que se presenta en [7].
swb2712: Un generador modificado de resta con préstamo, como se describe en [8]. Este generador es similar a un generador de Fibonacci retardado sumativo con retardos 27 y 12, pero está modificado para que tenga un periodo mucho más largo de aproximadamente $2^{1492}$ . El generador funciona de forma nativa con precisión doble para crear valores U(0, 1) y son posibles todos los valores en el intervalo abierto (0, 1). El algoritmo randn que se utiliza de forma predeterminada para las secuencias swb2712 es el algoritmo zigurat [7], pero con el generador swb2712 debajo.
Nota
Este generador es idéntico al que utiliza la función rand a partir de la Versión 5 de MATLAB, que se activa con rand('state',s).

Algoritmos de transformación

Inversion: Computa una variación aleatoria normal aplicando la función de distribución acumulativa inversa normal estándar a una variación aleatoria uniforme. Cada valor normal consume exactamente un valor uniforme.
Polar: El algoritmo polar de rechazo, como se describe en [1]. De media, cada valor normal consume aproximadamente 1,27 valores uniformes.
Ziggurat: El algoritmo zigurat, como se describe en [7]. De media, cada valor normal consume aproximadamente 2,02 valores uniformes.

Configuración de una secuencia

Una secuencia s de números aleatorios tiene propiedades que controlan su comportamiento. Para acceder o cambiar una propiedad, utilice la sintaxis p = s.Property y s.Property = p.

Por ejemplo, puede configurar el algoritmo de transformación para generar valores seudoaleatorios distribuidos normalmente cuando utiliza randn. Genere valores seudoaleatorios distribuidos normalmente con el algoritmo de transformación Ziggurat predeterminado.

s1 = RandStream('mt19937ar');
s1.NormalTransform

ans = 'Ziggurat'

r1 = randn(s1,1,10);

Configure la secuencia para utilizar el algoritmo de transformación Polar para generar valores seudoaleatorios distribuidos normalmente.

s1.NormalTransform = 'Polar'

s1 =
mt19937ar random stream
             Seed: 0
  NormalTransform: Polar

r2 = randn(s1,1,10);

Al generar números aleatorios con la distribución uniforme con rand, también puede configurar la secuencia para generar valores seudoaleatorios antitéticos, es decir, los valores normales restados de 1 para los valores uniformes.

Cree 6 números aleatorios con distribución uniforme a partir de la secuencia s.

s2 = RandStream('mt19937ar');
r1 = rand(s2,1,6)

r1 =
    0.8147    0.9058    0.1270    0.9134    0.6324    0.0975

Restablezca el estado inicial de la secuencia. Cree otros 6 números aleatorios con la propiedad Antithetic configurada como verdadera. Compruebe que estos 6 números aleatorios sean iguales a los números aleatorios generados de manera aleatoria restados de 1.

reset(s2)
s2.Antithetic = true;
r2 = rand(s2,1,6)

r2 =
    0.1853    0.0942    0.8730    0.0866    0.3676    0.9025

isequal(r1,1 - r2)

ans =
   1

En vez de configurar las propiedades de una secuencia una por una, puede guardar y restaurar todas las propiedades de una secuencia s utilizando A = get(s) y set(s,A), respectivamente. Por ejemplo, configure todas las propiedades de la secuencia s2 para que sean las mismas que en la secuencia s1.

A = get(s1)

A =
               Type: 'mt19937ar'
         NumStreams: 1
        StreamIndex: 1
          Substream: 1
               Seed: 0
              State: [625x1 uint32]
    NormalTransform: 'Polar'
         Antithetic: 0
      FullPrecision: 1

set(s2,A)

               Type: 'mt19937ar'
         NumStreams: 1
        StreamIndex: 1
          Substream: 1
               Seed: 0
              State: [625x1 uint32]
    NormalTransform: 'Polar'
         Antithetic: 0
      FullPrecision: 1

Las funciones get y set le permiten guardar y restaurar toda la configuración de una secuencia para que los resultados se puedan reproducir exactamente más adelante.

Restaurar estado de un generador de números aleatorios para reproducir la salida

La propiedad State representa el estado interno del generador de números aleatorios. Puede guardar el estado de la secuencia global en un punto determinado de su simulación al generar números aleatorios para reproducir los resultados más adelante.

Utilice RandStream.getGlobalStream para devolver un identificador a la secuencia global, es decir, la secuencia global actual de la que rand genera números aleatorios. Guardar el estado de la secuencia global.

globalStream = RandStream.getGlobalStream;
myState = globalStream.State;

Con myState, puede restaurar el estado de globalStream y reproducir resultados anteriores.

A = rand(1,100);
globalStream.State = myState;
B = rand(1,100);
isequal(A,B)

ans = logical
   1

rand, randi, randn y randperm acceden a la secuencia global. Dado que todas estas funciones acceden a la misma secuencia subyacente, llamar a una de ellas afecta a los valores que producirá el resto de ellas en las siguientes llamadas.

globalStream.State = myState;
A = rand(1,100);
globalStream.State = myState;
C = randi(100);
B = rand(1,100);
isequal(A,B)

ans = logical
   0

También puede restablecer una secuencia a su configuración inicial con la función reset.

reset(globalStream)
A = rand(1,100);
reset(globalStream)
B = rand(1,100);
isequal(A,B)

ans = logical
   1

Referencias

[1] Devroye, L. Non-Uniform Random Variate Generation, Springer-Verlag, 1986.

[2] L’Ecuyer, P. “Good Parameter Sets for Combined Multiple Recursive Random Number Generators”, Operations Research, 47(1): 159–164. 1999.

[3] L'Ecuyer, P. and S. Côté. “Implementing A Random Number Package with Splitting Facilities”, ACM Transactions on Mathematical Software, 17: 98–111. 1991.

[4] L'Ecuyer, P. and R. Simard. “TestU01: A C Library for Empirical Testing of Random Number Generators,” ACM Transactions on Mathematical Software, 33(4): Article 22. 2007.

[5] L'Ecuyer, P., R. Simard, E. J. Chen, and W. D. Kelton. “An Objected-Oriented Random-Number Package with Many Long Streams and Substreams.” Operations Research, 50(6):1073–1075. 2002.

[6] Marsaglia, G. “Random numbers for C: The END?” Usenet posting to sci.stat.math. 1999. Available online at https://groups.google.com/group/sci.crypt/browse_thread/ thread/ca8682a4658a124d/.

[7] Marsaglia G., and W. W. Tsang. “The ziggurat method for generating random variables.” Journal of Statistical Software, 5:1–7. 2000. Available online at https://www.jstatsoft.org/v05/i08.

[8] Marsaglia, G., and A. Zaman. “A new class of random number generators.” Annals of Applied Probability 1(3):462–480. 1991.

[9] Marsaglia, G., and W. W. Tsang. “A fast, easily implemented method for sampling from decreasing or symmetric unimodal density functions.” SIAM J. Sci. Stat. Comput. 5(2):349–359. 1984.

[10] Mascagni, M., and A. Srinivasan. “Parameterizing Parallel Multiplicative Lagged-Fibonacci Generators.” Parallel Computing, 30: 899–916. 2004.

[11] Matsumoto, M., and T. Nishimura.“Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudorandom Number Generator.” ACM Transactions on Modeling and Computer Simulation, 8(1):3–30. 1998.

[12] Matsumoto, M., and M. Saito.“A PRNG Specialized in Double Precision Floating Point Numbers Using an Affine Transition.” Monte Carlo and Quasi-Monte Carlo Methods 2008, 10.1007/978-3-642-04107-5_38. 2009.

[13] Moler, C.B. Numerical Computing with MATLAB. SIAM, 2004. Available online at https://www.mathworks.com/moler

[14] Park, S.K., and K.W. Miller. “Random Number Generators: Good Ones Are Hard to Find.” Communications of the ACM, 31(10):1192–1201. 1998.

[15] Salmon, J. K., M. A. Moraes, R. O. Dror, and D. E. Shaw. "Parallel Random Numbers: As Easy As 1, 2, 3." In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC11). New York, NY: ACM, 2011.