d. Generación de Muestras de Variables Aleatorias

250 Conceptos de Probabilidad, Variables Aleatorias
y Procesos Estocásticos en Redes de Comunicaciones
Marco Aurelio Alzate Monroy
Universidad Distrital F.J.C.
179
d. Generación de Muestras de Variables Aleatorias
Hemos insistido en que la aleatoriedad de un experimento se refiere a la imposibilidad de predecir
su resultado antes de ejecutar el experimento. La característica más sobresaliente de un programa de
computador es, precisamente, que los resultados son perfectamente predecibles, pues resultan de la
realización de un algoritmo específico ¿Cómo puede una máquina así generar muestras de variables
aleatorias? La verdad es que no lo hace, pues no puede lanzar dados o monedas. En vez de esto, se
usa un algoritmo determinístico para generar secuencias de números completamente determinísticas
que "se comportan" como si fueran aleatorias. Por ejemplo, considere los primeros 100000 dígitos
de la expansión decimal de  (3141592653 5897932384 6264338327 9502884197 1693993751
0582097494
5505822317
5665933446
4592307816
2535940812
1284756482
4062862089
8481117450
3378678316
9862803482
2841027019
5271201909
5342117067
3852110555
1456485669
9821480865
9644622948
2346034861
1328230664
9549303819
0454326648
7093844609
6442881097
2133936072
6024914127 3724587006 6063155881 7488152092 0962829254 etc.). La gráfica izquierda de la Figura
83 muestra cómo la distribución de cada dígito es perfectamente uniforme (cada dígito aparece
cerca de 10000 veces). Si condicionamos en el dígito anterior, por ejemplo si calculamos la
frecuencia relativa de los dígitos precedidos por un 0, como muestra la gráfica intermedia de la
Figura 83, la distribución sigue siendo uniforme (cada dígito aparece cerca de 1000 veces después
de un 0). Más aún, la función de autocovarianza normalizada es un impulso unitario, lo cual enfatiza
la sugerencia de independencia, pues ninguna sub-secuencia de dígitos parece contener información
sobre cuál será el siguiente dígito, como muestra la gráfica de la derecha de la Figura 83 . Otra
característica interesante es que una secuencia de N dígitos sólo se repetirá cerca de 105-N veces. Por
ejemplo, la secuencia '3' aparece 10026 veces, la secuencia '31' aparece 975 veces, la secuencia '314'
aparece 95 veces, la secuencia '3141' aparece 9 veces y la secuencia '31415' aparece 2 veces, pero la
secuencia '314159' sólo aparece una vez (de hecho aparecen todas las secuencias de 4 dígitos o
menos, pero sólo aparece el 37% de las secuencias de 5 dígitos). Si pudiésemos generar 100000
dígitos decimales realmente aleatorios, todas estas características también se cumplirían.
distribución de los dígitos de 
distribución de los dígitos de  que siguen a un 0
0.15
0.15
Covarianza entre los dígitos de 
1.2
0.1
0.05
0.8
0.1
Covarianza
fracción de ocurrencias
fracción de ocurrencias
1
0.6
0.4
0.05
0.2
0
0
0
1
2
3
4 5
dígito
6
7
8
9
0
0
1
2
3
4 5
dígito
6
7
8
9
-20
-10
0
separación
10
20
Figura 83. Características estadísticas de los primeros 100000 dígitos de la expansión decimal de 
Así como existen algoritmos para generar la secuencia de dígitos de , que es una secuencia
determinística que satisface todos los criterios de aleatoriedad, existen diferentes algoritmos para
generar secuencias de números enteros que satisfagan criterios de aleatoriedad como los
mencionados anteriormente, aunque en realidad sean completamente determinísticos. Por esa razón
se conocen como secuencias pseudo-aleatorias. La idea es que un generador de números enteros
entre 0 y m-1, que satisfaga criterios como los anteriores, puede usarse para generar números
250 Conceptos de Probabilidad, Variables Aleatorias
y Procesos Estocásticos en Redes de Comunicaciones
Marco Aurelio Alzate Monroy
Universidad Distrital F.J.C.
180
racionales en el intervalo [0,1) si dividimos cada entero de la secuencia entre m; estos racionales
podrían interpretarse como probabilidades en el programa simulador.
El área de la generación de números pseudo-aleatorios y la generación de muestras con otras
distribuciones distintas a la uniforme ha sido ampliamente estudiada desde la teoría de números y
ha arrojado muy interesantes algoritmos. Aquí mencionaremos unos pocos aspectos introductorios
para que el estudiante tenga una idea muy general y pueda ser, al menos, un buen usuario de las
funciones disponibles en las bibliotecas de los lenguajes de programación que utiliza.
Uno de los algoritmos más usados en la generación de secuencias uniformemente distribuidas de
números enteros es el generador congruencial lineal que, partiendo de una semilla Z0, genera los
términos de la secuencia así:
Z n   Z n 1    mod m , n  1
Este es un algoritmo rápido que requiere poca memoria. Como los números obtenidos se evalúan
módulo m, el conjunto de números que puede generar es Z{0,1,…,m-1}. Si la semilla Z0 se repite
después de K iteraciones, Z0 = ZK, toda la secuencia se repetirá a partir de ahí y, por consiguiente, el
período del generador será |Z| = K. La teoría de números describe condiciones en , , m y Z0 para
que (1) la secuencia esté uniformemente distribuida en Z, (2) los números de la secuencia sean
estadísticamente independientes, y (3) tengan un período tan grande como sea posible (ojalá m) –en
efecto, si un programa de simulación necesita 10000 números aleatorios, no sirve un generador que
repita la secuencia después de 5000 muestra–. Obsérvese que, si hace falta, cualquier secuencia
puede ser reproducible inicializando Z0 adecuadamente; esto es de gran importancia para comparar
resultados de diferentes modelos de simulación bajo idénticas condiciones.
Por ejemplo, empezando con Z0 = 0, el generador lineal congruencial {Zn = (5Zn-1 + 1) mod 16}
genera la secuencia {1, 6, 15, 12, 13, 2, 11, 8, 9, 14, 7, 4, 5, 10, 3, 0} que es perfectamente uniforme
y tiene máximo período, aunque se pueden ver fácilmente ciertas correlaciones (en cuatro ocasiones
ocurre el par (n, n+1), pero nunca ocurre el par (n, n+2), por ejemplo). Un generador útil debe tener
un valor alto de m (para un período suficientemente grande) y un valor alto de  (para proporcionar
suficiente aleatoriedad). Algunos generadores lineales congruenciales usados en diferentes
lenguajes de programación son los siguientes:
Zn = (1103515245 Zn-1 + 12345)
Zn = (22695477 Zn-1 + 1)
Zn = (6364136223846793005 Zn-1 + 1)
Zn = (16807 Zn-1 + 0)
etc.
mod 231 (la librería de GNU para el lenguaje C)
mod 232 (C/C++ de Borland)
mod 264 (la librería estándar C de newlib)
mod 231-1 (original de matlab®)
Sin embargo, existen otros algoritmos recientemente desarrollados con mejores características. Por
ejemplo, matlab permite escoger entre el generador multiplicativo congruencial anterior y otros
métodos como el siguiente generador multiplicativo retardado de Fibonacci,
250 Conceptos de Probabilidad, Variables Aleatorias
y Procesos Estocásticos en Redes de Comunicaciones
Marco Aurelio Alzate Monroy
Universidad Distrital F.J.C.
181
Zn = (Zn-31 * Zn-63) mod 264
que tiene un período de 2124, o el más aceptado en la actualidad, el "Mersene twister", que tiene un
increíble período de 219937-1 y mucho mejores características estadísticas, aunque requiere más
recursos de cómputo y de memoria para su evaluación. Para calcular Zn se usa una forma
generalizada de generadores por registros de desplazamiento y generadores retardados de
Fibonancci. Matlab usa, por defecto, el generador Mersene twister y le da al usuario la opción de
cambiarlo por otros de su preferencia.
Con cualquiera de estos generadores, la secuencia rn = Zn/m corresponde a una serie de números
uniformemente distribuidos en el intervalo [0, 1), independientes entre ellos, de manera que
podemos usar la secuencia como una generadora de probabilidades. En la mayoría de lenguajes de
programación, existe la función rand() para generar el siguiente número de la secuencia (como en
matlab). Así, por ejemplo, para simular la lanzada de una moneda, asignamos los valores 0 y 1 a los
elementos del espacio muestral cara y sello, respectivamente, y escogemos cada uno de los posibles
valores con probabilidad 1/2. Esto es, evaluamos r = rand() y, si r cae en el intervalo [0,0.5),
escogemos 0 y, si r cae en el intervalo [0.5,1), escogemos 1.
>> lado = {'cara','sello'};
>> moneda = (rand()>=0.5); disp(lado{moneda+1})
Cada vez que repitamos la última línea obtendremos el resultado 'cara' o 'sello' con idéntica
probabilidad, independientemente de anteriores resultados. Este algoritmo equivale a invertir la
función de distribución acumulativa de Bernoulli, como muestra la Figura 84.
FX(x)
1.0
r2 (sello)
 x  FX1  rnd () 
0.5
r1 (cara)
0
x
0 (cara)
1 (sello)
Figura 84. Obtención de muestras de variables de Bernoulli mediante inversión de la CDF
Si tenemos un enlace que permanece ocupado el 80% del tiempo y queremos generar una muestra
de cuántos enlaces están ocupados en un instante dado usamos la siguiente línea:
>> EnlacesOcupados = (rand()>=0.2)
que retorna 1 con probabilidad 0.8 y 0 con probabilidad 0.2. Si transmitimos un bit sobre un canal
con BER=10-5 y queremos generar una muestra de cuántos bits se dañaron durante la transmisión,
usamos la siguiente línea
250 Conceptos de Probabilidad, Variables Aleatorias
y Procesos Estocásticos en Redes de Comunicaciones
Marco Aurelio Alzate Monroy
Universidad Distrital F.J.C.
182
>> BitsErrados = (rand()<1e-5)
que retorna 1 con probabilidad 10-5 y 0 con probabilidad 1-10-5. En cada caso, hemos invertido la
CDF de la variable aleatoria de Bernoulli correspondiente. Si deseamos generar el estado de
ocupación de 10 enlaces o el estado de error de 100000 bits usaríamos líneas como las siguientes
>> EnlacesOcupados = (rand(10,1)>=0.2);
>> BitsErrados = (rand(100000,1)<1e-5);
Igual podemos hacer, por ejemplo, para lanzar cien dados:
>> r = rand(100,1); dado = (r>=0)+(r>=1/6)+(r>=1/3)+(r>=1/2)+(r>=2/3)+(r>=5/6);
Con la anterior instrucción, el vector dado será un arreglo de 100 números uniformemente
distribuidos en {1,2,3,4,5,6} e independientes entre ellos pues, efectivamente, hemos hecho la
inversión de la CDF, dado = F-1(rnd()), como muestra la Figura 85.
FX(x)
1
5/6
Primero se
genera una
probabilidad,
r0 = rnd ()
r0
2/3
1/2
1/3
1/6
x
1
2
3
4
5
6
x0  FX1  r0 
Después se invierte la CDF
Figura 85. Simulación de lanzar un dado mediante inversión de la CDF
Supongamos que queremos generar una muestra de una variable aleatoria geométricamente
distribuida en {0,1,2,…} con parámetro p, esto es, P[X=k] = pk(1-p). La CDF de esta variable es
n
FX (n)   p j (1  p)  1  p n 1 de manera que, si hacemos rand() = FX(n), obtenemos n = log(1j 0
rand())/log(p) – 1. Como muestra la Figura 86, debemos seleccionar el mínimo entero no menor a n
y, como 1-rand() es otro número probabilidad, podemos usar rand() en vez de 1-rand(). Por
ejemplo, las siguientes líneas en matlab generan 10000 muestras de una variable aleatoria
geométrica y comparan su histograma con la pmf P[X=k] = pk(1-p), como muestra la gráfica
izquierda de la Figura 87.
>> r = rand(10000,1);
250 Conceptos de Probabilidad, Variables Aleatorias
y Procesos Estocásticos en Redes de Comunicaciones
Marco Aurelio Alzate Monroy
Universidad Distrital F.J.C.
183
>> n = ceil(log(r)/log(0.8) - 1); mx = max(n);
>> h = hist(n,0:mx); h = h/sum(h);
>> stem(0:mx,h,'bo'); hold on; stem(0:mx,0.2*0.8.^(0:mx),'rx')
FX(x)
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
x
0
2
4
6
8
10
12
14
16
18
20
Figura 86. Generación de muestras de una variable aleatoria geométrica
Como último ejemplo, consideremos una variable aleatoria T exponencialmente distribuida con
promedio 1/, FT(t) = 1-e-t. Usando como probabilidad 1-rand() obtenemos t = -log(rand())/. La
siguiente secuencia de instrucciones de matalab® genera los tiempos de servicio de 10000 paquetes
distribuidos exponencialmente con promedio 2 ms y estima la pdf mediante el histograma, como
muestra la gráfica derecha de la Figura 87.
>>
>>
>>
>>
>>
>>
t = -log(rand(10000,1))/500;
[h,s] = hist(t,20);
d = s(2)-s(1);
f = h/sum(h)/d;
g = 500*exp(-500*s);
plot(s,f,'b-',s,g,'r--')
pmf de la variable geométrica
0.2
pdf de la variable exponencial
400
estimada
teórica
estimada
teórica
350
0.15
300
250
0.1
200
150
0.05
100
50
0
0
10
20
30
40
50
60
0
0
0.005
0.01
0.015
0.02
Figura 87. Distribuciones teóricas y empíricas, estimadas de 10000 muestras pseudo-aleatorias
Desafortunadamente no todas las funciones acumulativas de distribución son invertibles. Por
ejemplo, la Gaussiana ni siquiera tiene una expresión en forma cerrada. Entonces se usan otras
propiedades estadísticas como, por ejemplo, el hecho de que un punto al azar en un cuadrado
250 Conceptos de Probabilidad, Variables Aleatorias
y Procesos Estocásticos en Redes de Comunicaciones
Marco Aurelio Alzate Monroy
Universidad Distrital F.J.C.
184
unitario, (u1, u2), puede transformarse en dos muestras independientes de variables aleatorias
gaussianas con media cero y varianza 1 mediante la transformación
 x1 
cos(2 u2 ) 
2
log(u
)


1
x 
 sin(2 u ) 
 2

2 
que en matlab se consigue así (por ejemplo)
>> u = rand(2,1);
>> x1 = sqrt(-2*log(u(1)))*cos(2*pi*u(2));
>> x2 = sqrt(-2*log(u(1)))*sin(2*pi*u(2));
En efecto, despejando para u1 y u2 obtenemos
u1  e


 x12  x22 /2
u2 
1
x
tan 1  2 
 x1 
2
cuyo jacobiano es la es la función de densidad de probabilidad conjunta de dos variables normales
independientes,
u1
u1
x1
x2
1  (x12  x22 )/2
 1  x12 /2   1  x 22 /2 

 
e
e
e


u2
u2
2
 2
  2

x1
x2
Usando este método Box-Muller, como se le conoce, podemos verificar la función de densidad de
probabilidad de las muestras obtenidas (ver Figura 88):
>>
>>
>>
>>
>>
>>
>>
>>
>>
u = rand(10000,2);
x1 = sqrt(-2*log(u(:,1))).*cos(2*pi*u(:,2));
x2 = sqrt(-2*log(u(:,1))).*sin(2*pi*u(:,2));
[h,s] = hist([x1; x2],20);
d = s(2)-s(1);
f = h/sum(h)/d;
g = exp(-s.*s/2)/sqrt(2*pi);
plot(s,f,'b-',s,g,'r--')
legend('estimada','teórica')
0.5
pdf de la variable gaussiana N(0,1)
estimada
teórica
0.4
0.3
0.2
0.1
0
-4
-3
-2
-1
0
1
2
3
4
5
Figura 88. Estimación de la pdf N(0,1) a partir de muestras generadas mediante Box-Muller
250 Conceptos de Probabilidad, Variables Aleatorias
y Procesos Estocásticos en Redes de Comunicaciones
Marco Aurelio Alzate Monroy
Universidad Distrital F.J.C.
185
Como la CDF en un punto es el área debajo de la pdf desde - hasta ese punto, otro algoritmo típico
consiste en generar puntos uniformemente distribuidos en un rectángulo que contenga toda el área
debajo de la pdf. Si el punto cae debajo de la curva, su coordenada x es una muestra de la variable
aleatoria. Si no, se rechaza el punto seleccionado y se repite el procedimiento hasta obtener una
muestra válida. Este 'método de rechazo', como se le conoce, tiene algunas variantes interesantes
que buscan mayor eficiencia o mayor exactitud, como el método ziggurat.
En general, los métodos de generación de números pseudo-aleatorios y de muestras de variables
aleatorias con distintas distribuciones es un campo de la teoría de números que avanza
aceleradamente y que tiene importantes aplicaciones no sólo en simulación sino en criptografía, por
ejemplo (aunque en esa área existen requerimientos adicionales de seguridad). En este aparte sólo
queríamos mostrar una breve introducción a algunos de los métodos utilizados, en especial si el
estudiante requiere generar muestras de sus propias distribuciones empíricas o teóricas.
Matlab, en particular, no sólo permite seleccionar el algoritmo con que se obtienen muestras
uniformes en [0,1) sino que también permite escoger el algoritmo con que se obtienen muestras de
la variable normal a partir de las uniformes (ver en la documentación del matlab los detalles de la
instrucción RandStream()). Excepto por razones de compatibilidad con programas existentes, o por
la necesidad de usar múltiples flujos, es adecuado usar los algoritmos por defecto que matlab
propone (Mersenne twister para números uniformes en [0,1) y Ziggurat para números gaussianos
con media cero y norma 1) pues, hasta ahora, son los que mejores características presentan. De otro
lado, matlab incluye rutinas para generar muestras de una gran cantidad de distribuciones mediante
la instrucción random() –Beta, Binomial, 2, F, Exponencial, Gamma, Gaussiana, Geométrica,
hipergeométrica, lognormal, Pareto, Poisson, Rayleigh, t de student, uniforme continua o discreta,
Valor Extremo, Weibull, etc.–.
Por ejemplo, si en el programa simulador del Listado 7 cambiamos los tiempos entre llegadas
constantes por tiempos exponencialmente distribuidos con promedio 1/ y cambiamos los tiempos
de servicio constantes por tiempos exponencialmente distribuidos con promedio 1/, habremos
pasado de simular la cola D/D/n/k a simular la cola M/M/n/k. Aunque el cambio es mínimo en
cuanto a la modificación del listado del programa simulador, el cambio en el análisis de los
resultados de simulación es monumental, como veremos en el siguiente subtítulo.
e. Simulación de colas aleatorias: Análisis de resultados
Como mencionamos, para convertir el programa simulador de la cola D/D/n/k a un programa que
simule la cola M/M/n/k basta con cambiar la distribución de los tiempos entre llegadas y tiempos de
servicio, de determinística a exponencial. Las tres líneas relevantes del programa del Listado 7 son
las siguientes
31
37
59
proximaLlegada = reloj + 1/lambda;
proximaSalida(ss) = reloj + 1/mu;
proximaSalida(k) = reloj + 1/mu;
% Programa la próxima llegada
% Programa la próxima salida de este servidor
% Programa la próxima salida de este servidor
250 Conceptos de Probabilidad, Variables Aleatorias
y Procesos Estocásticos en Redes de Comunicaciones
Marco Aurelio Alzate Monroy
Universidad Distrital F.J.C.
186
que ahora deberán expresar lo siguiente:
31
37
59
proximaLlegada = reloj + random('exp',1/lambda);
proximaSalida(ss) = reloj + random('exp',1/mu);
proximaSalida(k) = reloj + random('exp',1/mu);
% Programa la próxima llegada
% Programa la próxima salida de este servidor
% Programa la próxima salida de este servidor
Igualmente valdría la pena cambiar ahora la primera línea para indicar que es la cola M/M y no la
cola D/D:
1 function [EN, EQ, ET, EW, PB, G] = ColaMMnk(lambda,mu,servidores,cupo,tiempoSimulacion)
Mientras que al simular la cola D/D/n/k basta con obtener una traza y tomar las estadísticas
temporales de esa traza, pues todas las simulaciones arrojarán la misma traza ya que no hay ningún
componente aleatorio en ellas, cada ejecución del programa simulador de la cola M/M/n/k arrojará
una traza diferente debido a la inclusión de muestras de variables aleatorias. Esto genera una
incertidumbre en los resultados de la simulación de la cola M/M/n/k que no existían con la
simulación de la cola D/D/n/k. Por ejemplo, si ejecutamos tres veces el programa simulador de la
cola D/D/1 con =0.8 y =1, obtendremos tres veces la traza observada en la tercera gráfica de la
Figura 78, pero si ejecutamos tres veces el programa simulador de la cola M/M/1 con los mismos
parámetros, obtendremos tres trazas diferentes, como muestra la Figura 89.
Número de paquetes en un sistema M/M/1 con =0.8 y =1
10
simulación 1
simulación 2
simulación 3
Número de paquetes
8
6
4
2
0
0
10
20
30
40
50
tiempo en segundos
60
70
80
90
100
Figura 89. Tres simulaciones distintas de un mismo sistema M/M/1
Este ejemplo ilustra muy bien el concepto de proceso estocástico: Cada ejecución del programa
simulador es un experimento aleatorio que genera como resultado una traza del número de paquetes
en función del tiempo. Con n simulaciones tenemos n trazas, {Ni(t), t0}, i=1,2,…n. Si hacemos n
simulaciones, podemos estimar el número promedio de paquetes en el sistema en un instante t0
cualquiera mediante el estimador
1 n
Eˆ [ N (t0 )]   N i (t0 )
n i 1
Si suponemos que el sistema alcanza la estacionariedad, podríamos considerar que la función de
valor medio llega a dejar de depender del tiempo:
n
1
Eˆ [ N ]  lim  N i (t0 )
t0  n
i 1
250 Co
onceptos de Probabilidad,
P
Variables
V
Aleeatorias
y Proceesos Estocástticos en Redess de Comuniccaciones
Marco Aureliio Alzate Moonroy
Universidadd Distrital F.JJ.C.
187
mar más y más
m muestras se va reduciiendo la incertidumbre dee la
Y, si la varianza ess finita, al tom
estimacción:
1 n
N i (t0 )

n  t0  n
i 1
E[ N ]  lim lim
Por otrro lado, con una sola trazza, las propieedades de erggodicidad y eestacionariedaad nos permitten
saber que
q basta con simular duraante un tiempo
o infinito paraa calcular
1 T
E  N   lim  N (t )dt
T  T 0
tal com
mo lo hacíamo
os con la colaa D/D/n/k. Claaro, con la coola D/D/n/k noo queremos eestimar la meddia
de un proceso esto
ocástico estaccionario y errgódico, com
mo con la coola M/M/n/k,, sino que sóólo
mos minimizaar el efecto deel período tran
nsiente.
querem
os veinte simu
ulaciones paraa generar veinnte trazas com
mo las de la F
Figura 89 {Ni(t),
Así pues, realicemo
t0}, i=
=1,2,…20, y calculemos el promedio deel número de paquetes en cada instante:
for i=
=1:20
% Co
orre 20 sim
mulaciones
[t
traza,EN,EQ,
,ET,EW,G,PB
B] = ColaMMn
nk(0.77,1,1
1,inf,1500);
; % genera una traza
dt
t = traza(2:
:end,1) - traza(1:endt
-1,1);
% ti
iempo entre eventos
ar
reas = traza
a(1:end-1,2
2).*dt;
% Ár
reas de los rectángulo
os
EN
N = cumsum(a
areas)./tra
aza(2:end,1)
);
% Pr
romedio acu
umulado hast
ta
pl
lot(traza(2:
:end,1),EN);
% ca
ada instant
te (Figura 80
0)
ho
old on
end
Si tomamos muestraas de EN cad
da segundo, po
odemos prom
mediar los veinnte promedioos temporales en
instanttes enteros para consideraar el promed
dio temporal (en la variaable tiempo) y el promeddio
probab
bilístico (en la variable deel espacio mu
uestral). El reesultado se aaprecia en la Figura 90 paara
=0.75
5 y =1.
Figuraa 90. Trazas obtenidas
o
de 20
2 simulacionnes y función de valor meddio
Como podemos nottar, cada trazaa (y el promeedio estadísticco de las trazzas) sufre un transiente anntes
bilidad, aunqu
ue el promedio estadísticoo se estabilizaa más rápido,, tendiendo m
más
de alcaanzar la estab
claram
mente a lo que parece ser ell verdadero vaalor esperadoo del número de paquetes een el sistema. El
tiempo
o de simulació
ón debe ser taan pequeño como
c
se puedda por razoness de eficienciia, pero tambiién
250 Conceptos de Probabilidad, Variables Aleatorias
y Procesos Estocásticos en Redes de Comunicaciones
Marco Aurelio Alzate Monroy
Universidad Distrital F.J.C.
188
debe ser tan grande como se pueda para asegurar que el sistema alcanza el estado estable en el que
queremos medir los promedios. De otro lado, a menos que empecemos a tomar estadísticas una vez
haya transcurrido el período transiente, debemos alargar el tiempo de simulación para que el efecto
de este período transiente sobre las estadísticas temporales se haga despreciable. En general, el
tiempo de simulación se escoge de acuerdo con algunas pruebas piloto como las de la Figura 90.
El otro aspecto importante en la planeación de experimentos consiste en la selección del número de
repeticiones que debemos hacer para asegurar cierto nivel de confianza en los resultados de la
simulación. Una alternativa es hacer una única simulación durante un tiempo increíblemente largo,
pues la ergodicidad nos asegura que así nos acercaremos indefectiblemente al verdadero valor
medio. El problema es que puede que el tiempo que se requiera sea excesivamente largo y, aun así,
no tendríamos manera de medir si ya estamos suficientemente cerca del verdadero promedio o si
debemos alargar el tiempo de simulación para acercarnos con un nivel de confianza suficiente.
Supongamos, entonces, que escogimos un tiempo total de simulación T y repetimos la simulación n
veces para obtener n trazas {Ni(t),0  t  T}, i=1,2,…n. Con cada simulación calculamos el
promedio temporal  Ni  1 T Ni (t )dt , i  1, 2,..., n  y finalmente calculamos el promedio estadístico
0

de los n resultados, Nˆ 

T
n
1
 Ni . Este valor será nuestro estimador del número promedio de
n i 1
paquetes en el sistema, pero ¿qué tan cerca o qué tan lejos podemos estar del verdadero valor
esperado?
Consideremos cada Ni (el promedio temporal de cada simulación) como una variable aleatoria cuyo
valor esperado es el verdadero promedio que queremos encontrar, m = E[N], y que tiene una
varianza finita, 2 (esta varianza será más pequeña entre mayor sea el tiempo de simulación, T, pero
hemos supuesto T fijo para todas las simulaciones). Sabemos entonces que el promedio Nˆ tiende a
tener una distribución normal con media m y varianza 2/n, y la variable normalizada Z = ( Nˆ m)/(/n) tiende a tener una distribución normal con media cero y varianza 1, N(0,1), a medida que
aumenta el número de simulaciones, n. El 100% percentil superior de Z es el valor z para el cual
P[Z>z] = , de manera que P[|Z|z/2]=1-, como muestra la Figura 91.
fZ ( z) 
1  z 2 /2
e
2
fZ ( z) 
1  z 2 /2
e
2
área  P   z /2  Z  z / 2   1  
área  P  Z  z   
z
z
área  P  Z  z / 2    / 2
área  P  Z   z /2    / 2
 z / 2
z /2
z
Figura 91. Conceptos de percentil (a) e intervalo de confianza (b)
En matlab® resulta muy fácil hallar los percentiles de la distribución normal. Por ejemplo, el valor
de z que hace P[Z>z]=0.05 (el 5% percentil) se calcula así:
250 Conceptos de Probabilidad, Variables Aleatorias
y Procesos Estocásticos en Redes de Comunicaciones
Marco Aurelio Alzate Monroy
Universidad Distrital F.J.C.
189
>> z_05 = norminv(1 - 0.05,0,1)
pues la instrucción norminv invierte la CDF de la distribución normal, de manera que
norminv(p,,) retorna el valor de x que hace FX(x) = P[Xx] = p si X~N(,).
Tomemos la última expresión y remplacemos Z por ( Nˆ - m)/(/n) para multiplicar por /n en los
tres términos dentro del argumento de la probabilidad:


z 
Nˆ  m
 z 
 z /2   P    /2  Nˆ  m   / 2   1  
P   z / 2  Z  z /2   P   z / 2 
/ n
n
n 



Restando Nˆ en los tres términos y multiplicando por -1, obtenemos la siguiente expresión. Ella
determina un intervalo alrededor de Nˆ donde se encuentra la verdadera media, m, con una
confianza de 1-:
z 
z 

P  Nˆ   /2  m  Nˆ   /2   1  
n
n 

El intervalo que delimita m, Nˆ z/2/n, se denomina "intervalo de confianza" y la cantidad 1- es
el "nivel de confianza", pues 1- es la probabilidad con que podemos asegurar que m se encuentra
en dicho intervalo. Por ejemplo, si queremos un nivel de confianza del 95% (=0.05), usamos el
percentil 2.5%, que es z0.025 = 1.96 (norminv(1-0.025,0,1)=1.96) y, con él, calculamos el
intervalo del 95% de confianza en términos de Nˆ , 2 y n.
Desafortunadamente la anterior expresión se basa en el conocimiento de la varianza, 2, cuando
apenas estamos tratando de estimar la media, m. Una alternativa es estimar 2 a partir de las
muestras de las simulaciones mediante la expresión
S2 
1 n
( N i Nˆ ) 2

n  1 i 1
pero, en este caso, la variable normalizada T = ( Nˆ - m)/(S/n) no es exactamente normal, sino que
tiene una distribución t de student con n-1 grados de libertad. La Figura 92 muestra la distribución t
de student con 1, 2, 5 e  grados de libertad, donde el último caso corresponde exactamente a la
variable N(0,1). Como era de esperarse, la distribución t es más dispersa que la distribución normal
debido a que el uso del estimador S2 en vez de la varianza 2 introduce mayor incertidumbre.
fT (t ; df )
0.4
1 grado de libertad
2 grados de libertad
5 grados de libertad
9 grados de libertad
0.35
0.3
 grados de libertad (N(0,1))
0.25
0.2
0.15
0.1
0.05
0
-4
-3
-2
-1
0
1
2
3
4
t
Figura 92. Comparación entre la distribución t de student y la distribución gaussiana
250 Conceptos de Probabilidad, Variables Aleatorias
y Procesos Estocásticos en Redes de Comunicaciones
Marco Aurelio Alzate Monroy
Universidad Distrital F.J.C.
190
Procediendo de la misma manera que cuando se conocía la varianza,
t
S
t
S

P  Nˆ   / 2;n 1  m  Nˆ   /2;n 1   1  
n
n 

donde t;df es el 100% percentil de la distribución t de student con df grados de libertad, esto es,
P[T>t;df ]= si T es una variable aleatoria con distribución t de student con df grados de libertad. El
percentil t;df es igualmente fácil de calcular con matlab:
>> t_alfa_df = tinv(1-alfa,df);
El programa del Listado 8 realiza tantas simulaciones como sean necesarias para satisfacer unos
requerimientos de confianza específicos. Por ejemplo, supongamos que queremos obtener un
intervalo del 95% de confianza para el número promedio de paquetes en una cola M/M/1 con
=0.75 y =1, que sea 5% del promedio estimado. La manera de invocar el programa sería así:
>> [m,v,delta]=SimulacionesMMnk(0.75,1,1,inf,8000,0.05,0.95)
Un posible resultado sería el siguiente, en el que se requirieron 18 simulaciones:
Con
Con
Con
Con
Con
Con
Con
Con
Con
Con
Con
Con
Con
Con
Con
Con
Con
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
simulaciones,
simulaciones,
simulaciones,
simulaciones,
simulaciones,
simulaciones,
simulaciones,
simulaciones,
simulaciones,
simulaciones,
simulaciones,
simulaciones,
simulaciones,
simulaciones,
simulaciones,
simulaciones,
simulaciones,
el
el
el
el
el
el
el
el
el
el
el
el
el
el
el
el
el
número
número
número
número
número
número
número
número
número
número
número
número
número
número
número
número
número
promedio
promedio
promedio
promedio
promedio
promedio
promedio
promedio
promedio
promedio
promedio
promedio
promedio
promedio
promedio
promedio
promedio
de
de
de
de
de
de
de
de
de
de
de
de
de
de
de
de
de
paquetes
paquetes
paquetes
paquetes
paquetes
paquetes
paquetes
paquetes
paquetes
paquetes
paquetes
paquetes
paquetes
paquetes
paquetes
paquetes
paquetes
es
es
es
es
es
es
es
es
es
es
es
es
es
es
es
es
es
2.8562
2.9041
2.9323
2.9419
2.9330
2.8847
2.8981
2.9213
2.9673
2.9848
2.9848
3.0277
3.0239
3.0022
3.0099
3.0000
2.9862
más
más
más
más
más
más
más
más
más
más
más
más
más
más
más
más
más
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
menos
menos
menos
menos
menos
menos
menos
menos
menos
menos
menos
menos
menos
menos
menos
menos
menos
243.646%
47.3829%
24.7332%
16.6842%
12.6756%
11.1464%
9.3492%
8.1855%
7.8915%
7.1108%
6.4123%
6.5325%
6.0101%
5.8067%
5.4117%
5.1205%
4.9235%
function [m,v,delta] = SimulacionesMMnk(lambda, mu, servidores, cupo, ...
tiempoSimulacion, intervalo, confianza)
% [m,v,delta]=N_Simulaciones_MMnk(0.75,1,1,inf,8000,0.05,0.95)
% Ejecuta no menos de 5 ni más de 50 simulaciones de una cola M/M/servidores/cupo
% con tasa de llegadas lambda y tasa de servicio mu. Cada simulación se ejecuta
% durante tiempoSimulación segundos. El número de simulaciones se detiene si los
% intervalos del 95% de confianza se alcanzan a estimar dentro de más o menos 5%
% del promedio estimado. Devuelve la media (m), la varianza (v) y el intervalo de
% E[N], E[Q], E[T], E[W], Gamma y Pb
%
suma1 = [0 0 0 0 0 0];
% Acumula las muestras de las 6 variables
suma2 = [0 0 0 0 0 0];
% Acumula los cuadrados de las muestras
numeroSimulaciones = 0;
% Cuenta las simulaciones que debe hacer
porcentaje = 100;
% Nivel de confianza
while ((porcentaje/100 > intervalo) * (numeroSimulaciones<50) + ...
250 Conceptos de Probabilidad, Variables Aleatorias
y Procesos Estocásticos en Redes de Comunicaciones
Marco Aurelio Alzate Monroy
Universidad Distrital F.J.C.
191
(numeroSimulaciones<2))
% no menos de dos ni más de 50 simulaciones, esperando satisfacer el
% intervalo de confianza solicitado
numeroSimulaciones = numeroSimulaciones + 1; % Una simulación más
[traza,EN,EQ,ET,EW,G,PB] = ColaMMnk(lambda, mu, servidores, cupo, ...
tiempoSimulacion);
suma1 = suma1 + [EN EQ ET EW G PB]; % Acumula las medidas
suma2 = suma2 + [EN^2 EQ^2 ET^2 EW^2 G^2 PB^2]; % y sus cuadrados
if numeroSimulaciones > 1
m = suma1/numeroSimulaciones; % Estima los promedios y las varianzas
v = (suma2 - suma1.*suma1/numeroSimulaciones)/(numeroSimulaciones-1);
percentil = tinv(0.5+confianza/2,numeroSimulaciones-1);
delta = percentil*sqrt(v/numeroSimulaciones);
porcentaje = 100*delta(1)/m(1);
disp(['Con ' num2str(numeroSimulaciones) ' simulaciones, el ' ...
'número promedio de paquetes es ' num2str(m(1)) ' más o ' ...
'menos ' num2str(porcentaje) '%']);
end
end
Listado 8. Este programa hace tantas simulaciones como sean necesarias para obtener
un intervalo de confianza y un nivel de confianza específicos .
Es interesante notar que un estudio de simulación que involucre aleatoriedad jamás podrá arrojar
resultados como "el número promedio de paquetes es 3", sino que sus resultados deben tener la
siguiente forma: "Según el estudio de simulación, con un nivel de confianza del 95% el número
promedio de paquetes se encuentra entre 2.9 y 3.1".

Download Report