[R-es] teoría de colas simulación

[Olivier Nuñez]

Una alternativa creo más eficiente:

MM1.4<-function(n=100,lambda=.8,mu=1) #n: número de clientes, lambda:  
intensidad de arribos, mu: intensidad de los tiempos de servicio
{
ltime = rexp(n,lambda) #tiempos entre dos llegadas exponencial
stime = rexp(n,mu) #tiempos de servicio exponencial
ltime[1]=0;stime[1]=0 #inicialización de la cola
lltime=cumsum(ltime) #instantes de llegadas al sistema
#instantes de salidas del sistema con politica FCFS (First Come First  
Serve)
sstime=rep(0,n)
for(t in 2:n) sstime[t]=max(sstime[(t-1)],lltime[t]) + stime[t]
##########
lmax=max(lltime)
time=c(lltime,sstime[(sstime<lmax)]) #instantes de saltos hasta la  
ultima llegada
events=data.frame(time=time,event=c(rep(1,n),rep(-1,(length(time)-n))))
Events=events[order(events$time),] #eventos ordenados
Ntime= cumsum(Events$event) #numero de clientes en el sistema en los  
instantes de saltos
out=cbind(tjump=Events$time[-1],size=Ntime[-1])
plot(out,type="l")
out
}

con

ptm <- proc.time();MM1.3(10000);proc.time()-ptm
    user  system elapsed
   1.647   0.439   2.099

mientras

ptm <- proc.time();MM1.4(10000);proc.time()-ptm
    user  system elapsed
   0.457   0.037   0.496
--  
____________________________________

Olivier G. Nuñez
Email: onunez en iberstat.es
Tel : +34 663 03 69 09
Web: http://www.iberstat.es

____________________________________




El 17/04/2010, a las 13:14, jose luis escribió:

> Buenas.
>
> Estoy intentando simular el modelo M/M/1 en R.  Y he encontrado el  
> siguiente código. Me gustaría saber si existe algún paquete que  
> pueda hacer esto, o si no, cómo vectorizar este código para que sea  
> más eficiente.
>
> Gracias
>
> MM1=function(n=100,lambda=.8,mu=1) #n: número de saltos, lambda:
> #intensidad de arribos, mu: intensidad de los tiempos de servicio
> {
> i=0
> tjump=rep(0,n)
> size=rep(0,n)
> size[1]=i
> for(k in 2:n)
>    {
>    if(i==0){mutemp=0}else{mutemp=mu}
>    time= rexp(1,lambda/(lambda+mutemp)) #tiempos entre saltos, con
> #distribución exponencial de parámetro (lambda+mutemp)
>    if(runif(1)<lambda/(lambda+mutemp)){i=i+1}else{i=i-1}
>    size[k]=i
>    tjump[k]=time
>    }
> tjump=cumsum(tjump) #tjump=tiempos acumulados de los saltos
> plot(tjump,size,pch=20,type="o")
> out=list(tjump=tjump,size=size) #size=tamaño del sistema
> out
> }
>
> _______________________________________________
> R-help-es mailing list
> R-help-es en r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es