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

[Carlos J. Gil Bellosta]

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
>