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
}
Hola, ¿qué tal?
Creo que esto está cerca de lo más eficientemente que puede
implementarse el código original:
MM1.2 <- function( n=100, lambda=.8, mu=1 )
#n: número de saltos, lambda:
#intensidad de arribos, mu: intensidad de los tiempos de servicio
{
quot <- lambda/(lambda + mu)
size <- rep(0,n)
tjump <- rep( quot, n )
jump <- sapply( runif( n ) < quot, function( x ) ifelse( x, 1, -1 )
)
for( k in 2:n )
{
if( size[k-1] == 0 ) tjump[k] <- jump[k] <- 1
size[k] <- size[k-1] + jump[k]
}
tjump <- cumsum( sapply( tjump, function( x ) rexp(1, x ) ) ) #
tiempos acumulados de los saltos
plot(tjump,size,pch=20,type="o")
list(tjump=tjump,size=size) #size=tamaño del sistema
}
Un saludo,
Carlos J. Gil Bellosta
http://www.datanalytics.com
On 04/17/2010 01:14 PM, jose luis wrote:> 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
>
Y esto podría ser mejor, incluso, en función de los parámetros de entrada:
MM1.3 <- function( n=100, lambda=.8, mu=1 )
#n: número de saltos, lambda:
#intensidad de arribos, mu: intensidad de los tiempos de servicio
{
quot <- lambda/(lambda + mu)
tjump <- rep( quot, n )
jump <- sapply( runif( n ) < quot, function( x ) ifelse( x, 1,
-1 ) )
jump[1] <- tjump[1] <- 1
foo <- function( x ){
tmp <- which( x == -1 )
if( length( tmp ) == 0 ) return( 0 )
tmp[1]
}
while( primer.menos.uno <- foo( cumsum( jump ) ) ){
jump[ primer.menos.uno - 1 ] <- tjump[ primer.menos.uno -1 ] <- 1
}
tjump <- cumsum( sapply( tjump, function( x ) rexp(1, x ) ) ) #
tiempos acumulados de los saltos
size <- cumsum( jump )
plot(tjump,size,pch=20,type="o")
list(tjump=tjump,size=size) #size=tamaño del sistema
}
(Especialmente, si la cola simulada es del Ayuntamiento de Madrid).
Carlos J. Gil Bellosta
http://www.datanalytics.com
On 04/17/2010 01:14 PM, jose luis wrote:> 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
>
Gracias Carlos. Eres un crack. Carlos J. Gil Bellosta escribió:> Y esto podría ser mejor, incluso, en función de los parámetros de > entrada: > > MM1.3 <- function( n=100, lambda=.8, mu=1 ) > #n: número de saltos, lambda: > #intensidad de arribos, mu: intensidad de los tiempos de servicio > { > > quot <- lambda/(lambda + mu) > > tjump <- rep( quot, n ) > > jump <- sapply( runif( n ) < quot, function( x ) ifelse( x, 1, > -1 ) ) > jump[1] <- tjump[1] <- 1 > > foo <- function( x ){ > tmp <- which( x == -1 ) > if( length( tmp ) == 0 ) return( 0 ) > tmp[1] > } > > while( primer.menos.uno <- foo( cumsum( jump ) ) ){ > jump[ primer.menos.uno - 1 ] <- tjump[ primer.menos.uno -1 ] <- 1 > } > > tjump <- cumsum( sapply( tjump, function( x ) rexp(1, x ) ) ) # > tiempos acumulados de los saltos > size <- cumsum( jump ) > > plot(tjump,size,pch=20,type="o") > list(tjump=tjump,size=size) #size=tamaño del sistema > } > > (Especialmente, si la cola simulada es del Ayuntamiento de Madrid). > > Carlos J. Gil Bellosta > http://www.datanalytics.com > > On 04/17/2010 01:14 PM, jose luis wrote: >> 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 >> > > _______________________________________________ > R-help-es mailing list > R-help-es en r-project.org > https://stat.ethz.ch/mailman/listinfo/r-help-es
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