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