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

[Carlos J. Gil Bellosta]

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
>