Versión: 26/11
. Al final se ven comentarios
sobre la sesión del V 26/11
LABORATORIO de PROBABILIDAD I , Curso 2004/5.
Código de las sesiones de proyección (C-15-102),
con comentarios y sugerencias.
Importante:
El código es un medio, como todo lenguaje; pensar sobre
todo en el significado de lo que hacemos.
Es bueno que el resultado sea sugerente y visualmente
atractivo; pero no caer en el culto del diseño.
Los ejemplos dados son ante todo provocaciones: examinarlos bien,
pero
¡tomar la iniciativa cuanto antes!
x = rand
(escribir eso en la CommandWindow y pulsar
'Enter')
El número que devuelve Matlab como valor de x podía haber
sido, "con igual IP,
cualquier 0<x<1"
Esa es la manera intuitiva de decir qué es la distribución UNIFORME(0,1),
pero da problemas:
- la definición rigurosa debe renunciar a hablar de "la IP
de cada valor x" que es necesariamente 0, y decir en su lugar:
"la IP de caer en dos intervalos de igual longitud
debe ser igual";
- lo que la función rand puede hacer
es algo más modesto: cumplir con esa definición hasta una
precisión razonable (que no se acerque al "grano" del
ordenador...); probemos si
lo hace bien; lo bueno del programa es que puede hacer muchas repeticiones
en un instante:
%--------------------------------- "distribución UNIFORME(0,1)"
x = rand(1000,1); plot(sort(x))
Vemos
esos 1000 números como ordenadas,
que
el sort ha ordenado
de menor a mayor.
¿Ha salido "lo que debe"? ¿Cómo "esperamos que estén
colocados"
1000 números con esta distribución? Respuesta: como
átomos en un cristal, de modo que en intervalos iguales haya "el
mismo número de ellos".
Los 1000 números unif que
producimos a continuación siguen al pie de la letra ese sistema:
i=1:1000; unif =(i-1/2)/1000; % ... y podemos
compararlos con los x :
plot(i, sort(x), i, unif, 'r')
Para
asegurarnos de haberlo entendido: si la IP de estar en el intervalo
(0,1/2) fuese 4 veces mayor que en el (1/2,1), pero UNIF en cada uno de
ellos, ¿cómo
habría que distribuir los "números modelo"?
El
plot
que resulte entonces tendrá diferente pendiente en esas dos
mitades, pero ¿en cuál de ellas más? Tratar de
razonar la respuesta, y luego traducirla a
código, a ver si todo cuadra.
La idea clave es: los valores son las ordenadas, y en las abscisas queremos frecuencia, así que
éstas deben estar "regularmente espaciadas", y son las ordenadas
las que
deberán ser más o menos densas...
Pero un momento, una distribución aproximadamente UNIF(0,1) se
conseguiría también con la sucesión 1/2,
1/4, 3/4,
1/8, 3/8, 5/8, ...
¿ POR QUE NO VALE ESO COMO 'rand' ??...
Claro, queremos también que los sucesivos valores
sean independientes ;
ese es quizá el concepto más profundo que tenemos entre
manos (con suerte, iremos
entendiéndolo a lo largo del curso); en el caso que nos ocupa,
podemos poner a prueba consecuencias
de esa idea;
por ejemplo, si seleccionamos valores de la lista x
con cualquier criterio, los
que les preceden
deberían seguir
teniendo la distribución del total:
%--------------------------------- consecuencias de la independencia:
i = find(x < x(1)); %
selecciona los índices de los x(i) que cumplen la
condición
% los anteriores
a esos en la lista, deben tener igual distribución que el total:
n=size(i); s=sort(x(i-1)); plot(i,s,i(n),s(n),'r')
% ... o que los anteriores a estos otros:
i = find(x > x(1)); n=size(i); s=sort(x(i-1));
plot(i,s,i(n),s(n),'r')
Pensar otros "tests de independencia" y ponerlos a prueba!
% Comentario sobre el
código:
¿por qué ' < x(1) ' en lugar de p.ej.'
< 0.5 ' ?
% porque en ese caso x(1) podría haber cumplido la
condición, ' 1 ' sería uno de los
índices ' i ',
% y entonces ' i-1' contendría al ' 0 ',
que no es un índice válido (habría que
suprimirlo...);
% por la misma razón usamos los anteriores: podemos
usar los siguientes
(es decir, i+1)
% si cambiamos la condición: x < x(1000) ...
Por suerte, el "generador seudo-aleatorio" rand (que también
tiene mi calculadora, y la mayoría ---buscar una tecla rand, o RAN
,...) nos
basta para generar "cualquier otra distribución", por ejemplo...
%--------------------------------- simular la MONEDA:
ca = (rand < 1/2) % el
paréntesis sobra, pero facilita la lectura!!
A una expresión así, que afirma algo cierto o falso (y evaluable por
Matlab), el programa le asigna un valor lógico =
1 ó 0, que puede usarse como un número ordinario para
cualquier operación posterior; por ejemplo, la linea que sigue
simula R=12 repeticiones del
experimento que consiste en "tirar N=20 veces una moneda equilibrada", produce
valores 0,1, y los ordena para dar una especie de "gráfico":
N=20; R=12; sort(rand(N,R)<1/2)
% Sugerencia: la
función ' bar ' es una buena
alternativa para hacer un gráfico muy parecido a esto.
%------------------ para números N, R mayores, mejor
usar gráficos:
N=80; R=100;
n = sum(rand(N,R)<1/2);
plot(n)
Este gráfico es la traducción obvia de lo anterior: las
ordenadas
n
sustituyen a las columnas de '1's.
Pero la información de en qué orden salieron los
distintos valores de
n = no. de caras, es
irrelevante (pensar por qué), de modo que mejor así:
figure; plot(sort(n)); grid
% se ha
usado ' figure ' para que no
desaparezca el gráfico anterior y poder compararlos;
% grid ha producido
la rejilla: se puede
apreciar cómo facilita el ver los datos.
Sin embargo, hay otro tipo de asuntos en los que el orden sí es relevante: por
ejemplo el de "cuán injusto puede ser el azar", manteniendo
largo tiempo a uno de "los dos
jugadores" por debajo de sus expectativas; para eso, tenemos que
"observar el desarrollo del juego"; en
consecuencia hemos llamado ' t ' a lo que es "el
paso del tiempo", o si se prefiere, el contador
de tiradas:
%--------------------------------- observar el "desarrollo
del juego":
p=1/2; N=80; t=1:N;
Empezamos
con la moneda equilibrada ( p=1/2 );
el código nos deja escoger
cada vez
como p la
ordenada del punto pinchado con el ratón, y ejecuta el bucle de
nuevo hasta que escojamos p<0:
while p>0
prop = cumsum(rand(1,N)<p)./t; % la
proporción de caras hasta la tirada t
plot(t,prop, t,p,'r', 0,[0,1])
[x,p] = ginput(1);
% escogemos el valor
de p para la vuelta
siguiente
end
La idea usada para
simular la moneda funciona igual para...
%---------------------------------
simular un DADO:
d=ceil(6*rand)
Armados
con esto, podemos por ejemplo atacar estos problemas de la Hoja 1:
%---------------------------------
problema H1.1
% IP de los valores 9 y 10 como suma de 3 dados
% ( las IP calculadas son:
IP(9) =25/216 =11.57 % , IP(10) =1/8 =12.5 % )
N=1e4; R=100;
% no.tiradas y no.repeticiones
for i=1:R
s=sum(ceil(6*rand(3,N)));
frec(i,:)=sum([s==9 ; s==10]')/N;
end
mf=mean(frec); frecuencias_medias = mf
plot(1:R,frec,
[0,R],[1;1]*mf,':k', 0,[0,0.15])
title('frecuencias de 9 y de 10
en 100 repeticiones de 10^4 tiradas')
xlabel(' en linea discontinua,
las frecuencias medias ')
%---------------------------------
problema H1.3.
% Se trata de un dado que se tira en total d veces, donde
% d= no. de puntos en la 1a.
tirada
% Queremos estimar la IP de B = "no.impar de seises en las d
tiradas"
R=1e4;
%
no.repeticiones de la simulación
for i=1:R
d=ceil(6*rand);
% 1a. tirada del dado
d(2:d)=ceil(6*rand(d-1,1)); % las
restantes d-1 tiradas
nim6(i)=mod(sum(d==6),2);
% indicatriz de B = "no.impar de seises"
end
format short;
p=1-mean(nim6) % la IP estimada de 'no B'
N=length(nim6);
remues=reshape(nim6,N/100,100);
sd_p=std(mean(remues))/10 % su SD
estimada por remuestreo
Se trata ahora de
flexibilizar un poco la idea:
%-------------------- simular una v.a. con unos pocos valores enteros
P=[.2 .3 .5]; cp=cumsum(P);
x = 1 + sum(rand > cp)
% toma valor k si p_1+...+p_k-1 < rand <
p_1+...+p_k
% luego lo toma con IP = p_k
for i=1:1e4
x(i) = 1 + sum(rand > cp);
end
frec=sum([x==1;x==2;x==3]')/1e4;
bar(1:3,frec); grid
Eso se puede usar
por ejemplo para el asunto de la...
%------------------------------ taquilla del teatro
% (ejemplo que se
propuso en el Lab el V 22/10)
P=[.3 .4 .1 .1 .05 .05];
% distribución del tamaño de un grupo
c=cumsum(P); Aforo=1000; R=100; B=30;
minb=1:R;
for i=1:R
n=0; b=B; mb=b;
while n < Aforo
k = 1 + sum(rand
> c); n=n+k;
if mod(k,2)==1, d=2*round(rand)-1;
% d = +-1 con igual IP
b=b+d; mb=min(mb,b); end
end
minb(i)=mb;
end
plot(sort(minb)); grid;
title('número
mínimo de billetes restantes (inicial=30)', 'Fontsize',14)
Sobre los asuntos
propuestos el V 5/11 :
Sobre el problema H1.12. (urna de
Polya):
Aunque la
pregunta de interés es "qué ocurre a la larga con la
proporción de azules",
no está de más observar un poco cómo ocurre, es decir,
observar la evolución de la
urna unas cuantas veces; para eso no es preciso hacer simulaciones muy largas, es mejor
ver cómo "se va estabilizando la composición" al comienzo
del proceso...
%------------- 12 simulaciones con urna
inicial [a0,b0]
T=100; a0=3; b0=2; %
duración y valores iniciales
t=1:T; n0=a0+b0;
for j=1:12
% no. de repeticiones
n=n0; a=a0;
% se inicializa la urna
for i=1:T
a=a+(rand<a/n); n=n+1;
% evolución, que se anota en la matriz 'at':
at(i,j)=a;
end
end
plot(n0+t,at,'k',t,t*a0/n0,':r'); % se
ha marcado en ':' la proporción inicial
axis( 'equal',[0,n0+T,0,a0+T])
title( 'composición de la urna en 8 repeticiones' ,
'Fontsize',16)
ylabel( ['número de azules; inicial = '
num2str(a0) ' de ' num2str(n0)] , 'Fontsize',14)
xlabel( 'número total' , 'Fontsize',14)
Pero para ver la distribución
de la proporción límite, conviene hacer más
repeticiones,
y más largas. No es muy difícil hallar fórmulas
para esa proporción límite, que se han
usado para compararlas con la distribución empírica:
%------------- 400 simulaciones, para ver
distribución del estado final
R=400; T=800; t=1:T;
for j=1:R
n=n0; a=a0;
for i=t, a=a+(rand<a/n); n=n+1; end
af(j)=a/n;
% se guarda la proporción final de cada
repetición
end
x=(1:R)/R; % 'x' es simplemente R puntos a
igual espacio sobre [0,1]
f=x.^(a0-1).*(1-x).^(b0-1);
% 'f' es la
función f(x) que da la densidad de IP para el lim(a/n),
% si la urna
tenía valores iniciales
a0, b0
F=cumsum(f)/sum(f);
% 'F' es
la correspondiente función de distribución;
% la cte de normalización que
falta en f se suple con la '/sum(f)'
figure; plot(sort(af),x,x,F,':'); axis square
xlabel('valor final de a/n en 400
repeticiones' , 'Fontsize',14)
Sobre el problema de la aguja de
Buffon:
No es
difícil simular el que la aguja corte a
las abscisas x=entero;
basta elegir la de un
extremo de la aguja como UNIFORME(0,1), el ángulo
como UNIFORME(0,2*pi), y
preguntar si el otro extremo está en (0,1), es decir, si tiene floor == 0 :
x=rand, t=2*pi*rand, nocorta=
floor(x+cos(t))==0
...o simular con la misma idea el que corte a
la retícula x=entero, y=entero :
x=rand, y=rand, t=2*pi*rand
nocorta= (floor(x+cos(t))==0)*(floor(y+sin(t))==0)
...pero la pregunta interesante
es:
¿ cuántas
repeticiones nos
harían falta para tener por ejemplo 3 dígitos fiables de
la IP ??
La respuesta es desalentadora: para que el "error
esperable" en la frecuencia observada
(su desviación típica) baje hasta < 0.5E-3 , hacen
falta, en números redondos, 10^6 reps,
y 100 veces más para cada "dígito fiable" adicional !!!
Por el contrario, hallar la IP(corte) =
2/pi es un sencillo ejercicio de Calculo II.
Veamos cómo sería el código para 4E6
repeticiones (que "garantizan" un error < 0.5E-3
con IP > 95 %) :
%------------ 4e6 tiradas de la aguja de Buffon
for i=1:2000
for j=1:2000
x=rand; t=2*pi*rand; nc(j)=(floor(x+cos(t))==0);
end
nocorta(i)=mean(nc);
end
IP_nocorta=mean(nocorta)
% pBuff estimada
Esto tendrá ocupado al ordenador un buen rato, y
debería dar un valor cercano a
pBuff=1-2/pi
% pBuff = 0.363380 (hallada
integrando)
pero el
tamaño esperable del error será la d.t. de la
frecuencia
observada:
sd=sqrt(pBuff*(1-pBuff))/2000
% = 2.4049e-004
de modo que sólo 3 cifras serán de fiar. Pero en
realidad hacer esto es una insensatez:
puesto que estamos simplemente sembrando el espacio muestral
[0,1] x [0,2*pi] con
4 millones de puntos, por qué no ponerlos en hileras ordenadas y dejar
dormir el 'rand' ??
(es
decir, hacer una
"integración numérica" mmmuy ingenua) :
%-------- versión NO-RAND del código
anterior
N=2000; x=(0.5:N)/N; x=(x'>0)*x;
nc=(floor(x+cos(2*pi*x'))==0);
IP_nocorta=mean(mean(nc))
err=IP_nocorta-pBuff
% err = -5.2276e-006 , casi 50 veces menor que la sd calculada antes !!!
% ... y
con mucho menos de 1/50
del esfuerzo de cálculo
En realidad sigue siendo una operación
inútil, puesto que el cálculo cerrado de la integral
es tan fácil... pero
si la hacemos, podríamos sacar algo más:
por ejemplo,
SUGERENCIA: hallar la
distribución de x,t
condicionada a que la aguja corte...
Con la misma
idea, veamos el caso de la cuadrícula
(cuidado, el peligro es salirse de las
dimensiones máximas permitidas por Matlab a una variable; por
eso ahora usamos un for) :
%--------- 64e6 tiradas de la aguja de Buffon sobre
cuadrícula
clear; N=400; x=(0.5:N)/N; x=(x'>0)*x;
for t=0.5:N
nc=(floor(x+cos(2*pi*t/N))==0).*(floor(x'+sin(2*pi*t/N))==0);
nocorta(t+0.5)=mean(mean(nc));
end
IP_nocorta2=mean(nocorta)
pBuff2=1-3/pi;
%
pBuff2 = 0.045070 (esto es lo que sale ahora integrando)
err=IP_nocorta2-pBuff2 % err
= 2.0409e-005
% ... pero si
hubieramos "simulado" (tardando ni se sabe), ...
sd=sqrt(pBuff2*(1-pBuff2))/8000
% sd = 2.5932e-005
Asuntos
presentados el V 26/11 :
En el código que se verá más abajo se explota una
de las ideas clave de la simulación probabilista: cómo
generar variables con cualquier distribución
1-dimensional a partir de la uniforme (de la 'rand', 'ALEATORIO()', etc.).
La idea es la siguiente: si la función de
distribución deseada
F(x) = P[ X <= x ]
tiene
inversa g , es decir g(F(x)) = x , llamemos X a la variable generada
del modo siguiente:
X = g(U) , donde U es Uniforme(0,1)
Como g
es creciente (igual que su inversa F), se tiene
P[ X<=x ] = P[ g(U)<=g(F(x)) ]
= P[ U<=F(x) ] = F(x)
,
por ser g creciente y U uniforme(0,1). Es decir, X tiene
función de distribución F(x).
Esa idea se ha usado en el problema del álbum de cromos
para producir distribuciones con "algunos cromos poco probables"; para
ello, antes de discretizar U
en la forma usual con ceil(n*U)
, se han "desuniformizado" los valores de U acercando su distribución a
la F(x)=x^2 ; para conseguirlo
se ha usado como g
una combinación lineal de las dos inversas: x=u, x=sqrt(u) .
El código que sigue
- produce 4 valores distintos de un parámetro a que controla la "falta de
uniformidad" (a=0 da la
uniforme), y genera en consecuencia cada vez la "sucesión de
números de cromo" que va a llegar al comprador;
- muestra su distribución empírica en un
gráfico, y anota en la CWindow
la "rareza" del cromo más raro: el cociente de su frecuencia
relativa por la esperada en el caso uniforme (que es 1/n, con n = no.de cromos = 50 en el
código);
- simula R=500 veces el
proceso de llenar el álbum, y anota el tiempo (es decir, el número
de cromos) consumido cada vez;
- tras repetir esto para los 4 valores de a , produce un gráfico con la
distribución de todos esos tiempos en los 4 casos, con las
medias empíricas de cada uno y la media esperada del caso
uniforme.
En vista de esas medias, se observa que el número total
de valores U usados
supera 4*R*700 =1.4e6 ;
obsérvese la idea utilizada para reducir la masa del
cálculo: como cada simulación individual usa en el peor
caso unos pocos miles de valores, se ha creado para empezar una lista 'rand' de longitud N = n*R , que se "rota" tras cada
repetición para usarla circularmente como "universo muestral"
completo.
Hay extensiones obvias de este código en la
dirección de "cambiar cromos" y evaluar el ahorro que se
consigue con ello, pero dejamos esto a la iniciativa de quien quiera
hacerlo.
%-----------------------
simulación repetida del llenado de un álbum,
%----------------------- con 4
distribuciones cada vez menos uniformes
n=50; % no.de cromos
R=500; % no.de repeticiones
%-----------------------
T_esperado_uniforme =
n*sum(1./(1:n)) % = 224.96
si n = 50
N=n*R;
% N = tamaño del
universo muestral;
p=(1:N)/N;
% ordenadas para el subplot(2,1,1)
u=rand(N,1);
s=sqrt(u); % u es uniforme; s tiene
distribución F_s = x^2
for j=1:4
a=(j-1)/4;
% parámetro de no-uniformidad
x=a*s+(1-a)*u;
% la distribución tiene F'(0)=0 si a>0
k=ceil(n*x);
% lista
"circular" de N cromos
rareza_del_1 =
n*sum(k==1)/N % muestra la "rareza" de k=1
subplot(2,1,1);
plot(sort(k),p); hold on
for
i=1:R
% repeticiones de la
simulación
t=0;
f=ones(1,n); % tiempo= 0, faltan todos
los cromos
while
sum(f)>0 %
mientras falte alguno...
t=t+1; f(k(t))=0;
end
T(i,j)=t;
% anotar el tiempo T
empleado esta vez,
k=k([t+1:N ,
1:t]); % y "rotar" la lista
end
end
axis([0,n+1,0,1]); hold off;
% cierra el gráfico de distribuciones
title(['F de
distribución para los '
num2str(n) ' cromos'],'Fontsize',14)
subplot(2,1,2);
T=sort(T); Tm=mean(T);
fr=(1:R)/R;
plot(T,fr,'k',Tm,1/2,'r+',T_esperado_uniforme,1/2,'og');
axis([0,max(T(:))+n,0,1]);
title(['T de espera en ' num2str(R) ' repeticiones
(las "+" son valores medios)'],'Fontsize',14)
xlabel('cromos consumidos','Fontsize',14)
%=================== OUTPUT en CW, y
otros valores (en una ejecución):
% T_esperado_uniforme = 224.9603
% rareza_del_1 = 0.9900
% rareza_del_1 = 0.2160
% rareza_del_1 = 0.0720
% rareza_del_1 = 0.0280
% valores de T:
% min
=
119
116
207 600
% mediana
= 217
299
482 1682
% max(T)
= 437
1092 2636 3470
% Tm = 1000 * [ 0.2235 0.3336
0.7359 1.7859 ]
%==============================
En el siguiente ejemplo se ha vuelto a usar la idea de
generar X usando
la inversa de F(x) ; se
trata ahora de una cola en una ventanilla, y se han generado, para un
número prefijado de clientes, dos distribuciones:
- la de los intervalos entre sus tiempos de llegada, a los que
se da distribución
F(x) = 1 - exp(-cx)
usando la inversa:
X = g(U) = - log(1-U)/c
- la de sus "tiempos de servicio" en ventanilla, generados como
S
= exp(acos(1-2*U))
para dar una densidad en forma de "tarta estirada por la derecha", que
es frecuente en este tipo de datos (para servicios "rutinarios" como
compras de entradas, check-in en un aeropuerto, ...); se reescalan para
darles media prefijada s = 3 min
, y se muestran en un histograma.
La primera distribución es la Exponencial de parámetro c,
cuyas sumas parciales (los tiempos de llegada de los clientes)
serán lo que se llama un proceso
de Poisson de intensidad c , uno de los modelos clave del
Cálculo de Probabilidades; el parámetro c resulta ser el número
esperado de llegadas/1dad de tiempo, por lo que se ha escogido c = 0.3/min (si escogemos s*c > 1 , el servicio
estará sobresaturado y es de esperar que la cola crezca!!)
El resto del código
- produce el proceso mismo, anotando los tiempos A de acceso a ventanilla y B de fin de servicio, así
como el número de orden n
que cada cliente tiene en la cola al ingresar (=0 si encuentra la
ventanilla libre);
- calcula estadísticas y produce un gráfico con
los tiempos de llegada y acceso de cada cliente; puede observarse en
él un caso del principio general de "contar de dos maneras": la suma de
los minutos de cola esperados por
cada cliente es también la suma de las longitudes de cola observadas en cada
minuto (ambas miden "el área entre las dos curvas del
gráfico", en unidades "cliente*min");
eso permite usar la media de los tiempos de espera A-T para calcular la media de dichas
longitudes, que es diferente de la media de las longitudes n observadas a la llegada de
cada cliente, y de la mediana de esas mismas, también calculada:
%-----------------------
simulación de una cola de ventanilla
% ---------------------- N clientes;
% sus llegadas son un proceso de Poisson de intensidad c (en
unidades 1/min )
% tiempo medio de
servicio, s (en min )
N=90;
%---------------- producir intervalos
entre llegadas (en min)
c=0.3;
% para tener F(x) = 1 - exp(-cx) ,
tomamos:
U=rand(1,N); X = - log(1-U)/c;
%---------------- ver la
distribución producida, comparada con la teórica
subplot(2,1,1);
sx=sort(X); F= 1 - exp(-c*sx);
p=(1:N)/N; % p = ordenadas en [0,1]
plot(sx,p,sx,F,'r:');
axis([0,log(N)/c,0,1])
xlabel('minutos','Fontsize',14)
title(['Intervalos entre
clientes: Exp_c , con c = '
num2str(c)],'Fontsize',14)
%---------------- tiempos de servicio
(en min), con media:
s=3;
U=rand(1,N); S=exp(acos(1-2*U));
S=s*S/mean(S);
subplot(2,1,2); hist(S,30);
xlabel('tiempos de servicio, en
minutos','Fontsize',14)
%-------- PROCESO:
% tiempos T = llegada, A = acceso, B = fin servicio;
n= longitud cola
T=cumsum(X); A=T(1);
for i=1:N-1
n(i)=sum(A>T(i));
B(i)=A(i)+S(i); A(i+1)=max(T(i+1),B(i));
end
%-------- estadísticas
T_total = A(N)+S(N)
t_medio_en_cola = mean(A-T)
long_media_cola_por_min =
sum(A-T)/T_total
n=sort(n);
mediana_cola = n(N/2),
long_media_cola_por_clientes = mean(n)
%-------- gráfico
esc=[1;1]*(1:N); esc=[0,esc(1:2*N-1)];
% ordenadas para "escalera"
t=[T;T]; a=[A;A];
%
abscisas: tiempos de llegada y acceso
subplot(2,1,2);
plot(t(:),esc,'r', a(:),esc,'g');
axis([0,T_total,0,N])
xlabel('minutos','Fontsize',14); ylabel('clientes','Fontsize',14)
legend('Linea roja = tiempos de
llegada' , 'Linea verde =
acceso a ventanilla')
%=================== OUTPUT en CW (en
una ejecución):
% T_total = 313.5652
% t_medio_en_cola = 8.2583
% long_media_cola_por_min = 2.3703
% mediana_cola = 2
% long_media_cola_por_clientes =
2.9101
%
Al contrario del ejemplo anterior, que es diabólico
tratar de imitar con Excel, en
este caso es muy fácil replicar el código paso por paso,
como puede verse en este ejemplo.
También en este caso se ha iniciado una de las variantes
obvias: que haya v>1
ventanillas, con el sistema civilizado de que cada cliente guarda su
turno de llegada y ocupa la primera que quede libre; nótese que
bastan ligeros cambios de código, y que desde luego no hace
falta individualizar las diferentes ventanillas; no se han seguido
otras lineas obvias, como hacer estadísticas de tiempo ocioso en
ventanilla, tamaño de los intervalos ociosos,...
Se han añadido también m clientes que esperaban ya al
abrir la oficina:
%-------- VARIANTES: cola inicial, y
más ventanillas
m=8; % no. de
madrugadores
v=2; % no. de
ventanillas
k=1:v; Sv = S*v; X(1:m)=0*(1:m);
%-------- tiempos de llegada y de
acceso, y longitud cola
clear A; T=cumsum(X);
A(k)=T(k); B(k)=A(k)+Sv(k);
for i=v:N-1
n(i)=sum(A>T(i));
B(i)=A(i)+Sv(i);
A(i+1)=max(T(i+1),min(B(i+1-k)));
end
%-------- estadísticas
T_total = A(N)+S(N)
t_medio_en_cola = mean(A-T)
long_media_cola_por_min =
sum(A-T)/T_total
n=sort(n);
mediana_cola = n(N/2),
long_media_cola_por_clientes = mean(n)
%-------- gráfico
t=[T;T]; a=[A;A];
% abscisas: tiempos de llegada y acceso
figure
plot(t(:),esc,'r', a(:),esc,'g');
axis([0,T_total,0,N])
title(['Evolución con ' num2str(v) ' ventanillas y ' num2str(m) ' madrugadores'],'Fontsize',14)
xlabel('minutos','Fontsize',14)
legend('Linea roja = tiempos de
llegada' , 'Linea verde =
acceso a ventanilla')
%=================== OUTPUT en CW (en
una ejecución)
% con 2
ventanillas, y el tiempo de servicio duplicado:
% T_total = 292.7892
% t_medio_en_cola = 4.8466
% long_media_cola_por_min = 1.4898
% mediana_cola = 1
% long_media_cola_por_clientes =
1.8652