sábado, 27 de noviembre de 2010

Sistemas de colas en R. Programación

####Sistemas de colas
 
#1)Sistema de servicio discreto
#se considera un tiempo discreto
#se considera un servidor de una sola cola
#la llegada de los clientes al sistema se realiza según una secuencia independiente {a_i, i=1,...,k} donde a_k es el número de llegadas que se han producido hasta el instante k
#los clientes esperan hasta que el servicio queda libre
#el número de clientes en el sistema al tiempo k, está dado por n(k)=n(k-1)+a(k)-I{n(k-2)+a(k)>=1}, k>=2 (condición inicial: n(1)=0)
#esta recursión define la cadena de Markov {n(k)}, k>=1
#En matlab:
# m = 200; p = 0;2;
# N = zeros(1,m);
# A = geornd(1-p,1,m);
# for n = 2 : m
# N(n) = N(n-1) + A(n) - (N(n-1) + A(n) >= 1);
# end
# stairs((0 : m-1);N);
 
colas_discreto<-function(m=100,p=.2)
{
N=rep(0,m)
A=rgeom(m,1-p) #Aplica el modelo de secuencia de llegadas (En este caso geométrico: genera m números aleatorios geométricos, con probabilidad 1-p)
for(n in 2:m)
{
N[n]=N[n-1]+A[n]-(N[n-1]+A[n]>=1) #I{n(k-2)+a(k)>=1} es la función indicadora de distribución cero-uno
}
x=0:(m-1)
par(mfrow=c(1,2))
plot(N~x,pch=20,type="o")
plot(cumsum(N),pch=20,type="b")
out=list(N=N,cumsumN=cumsum(N))
out
}
 
#2) modelo M/M/1
#se considera un tiempo continuo
#se considera un servidor de una sola cola
#la llegada de los clientes al sistema se realiza según un proceso de Poisson con intesndidad lambda.
#los tiempos de servicio presentan distribución Poisson de parámetro mu. Por tanto, el servidor necesita un tiempo aleatorio distribuido exponencialmente para cada cliente, con una media de tiempo de servicio de 1/mu.
#el resultante tamaño del sistema N(t), t>=0, es un proceso de Markov en tiempo continuo
#in Matlab:
#if i==0 mutemp=0; else mutemp=mu; end
#The main loop just checks for jump up or jump down:
#if rand<=lambda/(lambda+mutemp)
# i=i+1; %jump up: a customer arrives
#else
# i=i-1; %jump down: a customer is departing
#end %if
#x(k)=i; %system size at time i
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=-log(runif(1))/(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
}
 
 
#modelo M/D/1
#el sistema presenta arribos con distribución Poisson de intensidad lambda(al igual que el M/M/1) pero cada tiempo de servicio se fija con largo=1, por lo que no es un Proceso de MArkov
#las partidas de los clientes ocurre en los tiempos: u_k=k+max(t_1,t_2-2+1,...,t_k-k+1) donde t_1,t_2,...,t_k son los tiempos de arribo del proceso de Poisson
#el tamaño del sistema N(t), t>=0, presenta saltos hacia arriba en cada t_k y hacia abajo en cada u_k.
#En matlab:
# arrsubtr=arrtime-(0:n-1)';
# arrmatrix=arrsubtr*ones(1,n);
# deptime=(1:n)+max(triu(arrmatrix))
# B = [ones(n; 1) arrtime; ¡ones(n; 1) deptime´];
# Bsort = sortrows(B, 2);
# order
# jumps = Bsort(:, 1); jumptimes = [0; Bsort(:, 2)];
# X = [0;cumsum(jumps)];
#library(gdata)
MD1=function(lambda=.95,tmax=100) #tmax=intervalo de simulación;lambda=intensidad de arribos
{
arrtimes=-log(runif(1))/lambda #arribos con distribución Poisson (tiempos entre arribos con ditribución exponencial)
i=1;while(min(arrtimes[i])<=tmax){arrtimes=c(arrtimes,arrtimes[i]-log(runif(1))/lambda);i=i+1}
n=length(arrtimes)
arrsubtr=arrtimes-0:(n-1) #t_k-(k-1)
arrmatrix=arrsubtr%*%t(rep(1,n))
deptime=(1:n)+max(upperTriangle(arrmatrix)) #max(arrtimes[1:n])??? #tiempos de salidas de los clientes u_k=k+max(t_1,...t_k-k+1)
col1=c(rep(1,n),rep(-1,n)); col2=c(arrtimes,deptime); dd=data.frame(col1,col2)
Bsort=with(dd,dd[order(dd[,2]),]) #saltos ordenados
 
jumps=Bsort[,1]
jumpstimes=c(0,Bsort[,2]) #tiempos de arribos o salidas
size=c(0,cumsum(jumps)) #tamaño del sistema
time=deptime-arrtimes #tiempos del sistema
par(mfrow=c(1,2))
plot(time,type="o",pch=20)
plot(size,type="o",pch=20)
out=list(jumps=jumps,jumpstimes=jumpstimes,size=size,time=time)
out
}
 
#3) modelo M/G/1 (generalización del M/D/1 al caso de tiempo de servicio con distribución general de media 1/mu)
#los arribos presentan distribución Poisson con intensidad lambda
#los tiempos de servicio son uniformes
#es una generalización del modelo M/D/1, para un tiempo de servicio con distribución general S y media 1/mu.
#como caso particular: cuando los tiempos de salida se distribuyen según una exponencial de parámetro mu, se tiene el modelo M/M/1
 
MG1=function(tmax=100,lambda=.9) #tmax: intervalo de simulación, lambda: intensidad de arribos
{
arrtime=-log(runif(1))/lambda #arribos con distribución Poisson
i=1
while(min(arrtime[i])<=tmax)
{
arrtime=c(arrtime,arrtime[i]-(log(runif(1))/lambda))
i=i+1
}
n=length(arrtime) #tiempos entre arribos
time=2*runif(n) #tiempos de servicio s_1,...,s_k
cumsumtime=cumsum(time)
arrsubtr=arrtime-c(0,cumsumtime[1:(n-1)]) #t_k-(k-1)
arrmatrix=arrsubtr*matrix(1,1,n)
deptime=cumsumtime+max(arrmatrix) #tiempos de salida u_k=k+max(t_1,...,t_k-k+1)
col1=c(rep(1,n),rep(-1,n)); col2=c(arrtime,deptime); dd=data.frame(col1,col2)
Bsort=with(dd,dd[order(dd[,2]),]) #saltos ordenados
jumps=Bsort[,1]
jumpstimes=c(0,Bsort[,2]) #tiempos de arribos o salidas
size=c(0,cumsum(jumps)) #tamaño del sistema
time=deptime-arrtime #tiempos del sistema
par(mfrow=c(1,2))
plot(time,type="o",pch=20)
plot(jumpstimes,size,type="o",pch=20)
out=list(jumps=jumps,jumpstimes=jumpstimes,size=size,time=time)
out
}
#4)modelo M/G/infinito (cada cliente tieme acceso a su propio servicio, no hay clientes en cola, solo se simulan los tiempos de llegada y servicios)
#En matlab: cuyos tiempos de servicio se distribuyen según una distribución de Pareto
# alpha = 1,5;
# servtimes = rand^(-1/(alpha-1))-1;
# servtimes = [servtimes; rand(npoints-1,1).^(-1/alpha)-1];
 
ParetoST=function(alpha=1.5,npoints=100) #generar los tiempos de servicio con distribución Pareto
{
servtimes=runif(1)^(-1/(alpha-1))-1 #generar una muestra uniforme y aplicar el cdf inverso (cdf: F(x) = 1-(1+x)^(-alpha) -Pareto-)
servtimes=c(servtimes,runif(npoints-1)^(-1/alpha)-1)
plot(servtimes,pch=20,type="o")
servtimes
}
 
 
#los arribos son procesos de Poisson homogéneos con intensidad lambda
#los tiempos de servicio tienen distribución de Pareto
#cada cliente tiene su propio servidor (no es un sistema de cola -queue-)
MGinfinito=function(tmax=100,lambda=1,alpha=1.5,npoints=100) #tmax: intervalo de simulación, lambda: intensidad de arribo
{
arrtimes=-log(runif(1))/lambda #arribos con distribución Poisson (tiempos entre arribos con ditribución exponencial)
i=1;while(min(arrtimes[i])<=tmax){arrtimes=c(arrtimes,arrtimes[i]-log(runif(1))/lambda);i=i+1}
npoints=length(arrtimes)
if(npoints>0){arrtimes=sort(runif(npoints)*tmax)} else {arrt=NULL}
servt=runif(1)^(-1/(alpha-1))-1 #Tiempo de servicio de Pareto,proceso de renovación estacionario
servt=c(servt,runif(npoints-1)^(-1/alpha)-1)
servt=10*servt #se elige la media de manera arbitraria
dept=arrtimes+servt #tiempos de salida
 
col1=c(rep(1,npoints),rep(-1,npoints)); col2=c(arrtimes,dept); dd=data.frame(col1,col2)
Bsort=with(dd,dd[order(dd[,2]),]) #saltos ordenados
jumps=Bsort[,1]
jumpstimes=c(0,Bsort[,2]) #tiempos de arribos o salidas
size=c(0,cumsum(jumps)) #tamaño del sistema
time=servt #tiempos del sistema?
 
par(mfrow=c(1,2))
plot(jumpstimes,type="o",pch=20)
plot(jumpstimes,size,type="o",pch=20)
out=list(jumps=jumps,jumpstimes=jumpstimes,size=size,times=servt)
out #jumptimes:tiempos de cambios de estados en el sistema,size: número de clientes en el sistema
}
 
 
#####Ejercicio
#Generar el tamaño de un sistema en el caso de un M/G/infinito estacionario en el intervalor [0, 5),
#con intensidad de llegada lambda=2 y siendo los tiempos de servicio variables de Paretto normalizadas con parámetros alpha=1.6, gamma=2.
#Versión estacionaria del M/G/infinito, con tiempos de servicio arbitrarios
#En el tiempo 0 existe:
# un número de clientes según la distribución Poisson(lambda*mu)
# tiempos de servicio según una distribución estacionaria (en nuestro caso distribución de Pareto normalizada)
#la versióne stacionaria del M/G/infinito, se logra mediante:
# generar las llegadas "negativas" aún activas en el tiempo 0 y los tiempos de servicio estacionarios.
# distrstat.m para generar un número aleatorio desde la distrib estacionaria
# tiempo=0, número de clientes=Poisson, tiempo de servicio=estacionario: nstart=rpois(1,lambda*servmu)
# si nstart>0, obtener la distribución y parámetros estacionarios statdist/statpar, generar los parámetros para el generador de números aleatorios rndpar1, obtener los tiempos de servicio estacionarios stattimes, adjudicar tiempo de arribo 0 a los clientes que llegan antes de t=0 arrtimes
# una vez que creamos el proceso de conteo, eliminaremos los 0 extra en el comienzo de los arribos "negativos", y agregaremos el tmax a los saltos para obtener gráficos correctos
 
#genera una matriz de númneros aleatorios según la distribución de Pareto normalizada
# pdf f(x) = alpha*gamma/(1+gamma*x)^(1+alpha), x>0,
# cdf F(x) = 1-(1+gamma*x)^(-alpha)
ParetoNorm=function(M=1,N=1,alpha=1.6,gamma=2) #M y N son filas y columnas de la matriz a crear, aplha es el parámetro de cola en la distribución, y gamma el parámetro de la distribución (permite el control del valor esperado. Para alpha>1 el valor esperado existe y es igual a 1/(gamma*(alpha-1)))
{
sample=((1-matrix(runif(M*N),M, N))^((-1/alpha)-1)/gamma)
sample #matriz MxN de números aleatorios
}
 
####chequear el script de MGinfinito?????!!!!
MGinfinito_Estac=function(tmax=5,lambda=2,alpha=1.6,gamma=2) #tmax: intervalo de simulación, lambda: intensidad de arribos, alpha y gammma son los parámetros para la distrib ParetoNorm
{
servmu=1
#servpar=c(alpha,1/(servmu*(alpha-1))) #parámetros de la distribución de tiempos de servicio #en nuestro caso tenemos "gamma"
servdist=ParetoNorm(M=1,N=1,alpha=alpha,gamma=gamma) #distribución de los tiempos de servicio (función externa para la generación de números aleatorios a partir de la distribución deseada, en nuestro caso: ParetoNorm)
servmu = 1/(gamma*(alpha-1)) #switch func2str(distr)?? #???"dpar" #extrae los valores esperados de los tiempos de servicio; caso "ParetoNorm(alpha,gamma)"
nstart=-log(runif(1))/(lambda*tmax);i=1;while(min(nstart[i])<=tmax){nstart=c(nstart,nstart[i]-log(runif(1))/(lambda*tmax));i=i+1}
npoints=length(nstart)
#se asume que el sistema es estacionario desde el tiempo "-infinito". En el tiempo 0 el número de clientes es según un número de Poisson, con tiempos de servicios estacionarios
if(nstart>0)
{
statdist=servdist #switch func2str(distr)?? #distribución estacionaria G(t)=1/mu*int_0^t(1-F(s))ds, para la distribución "ParetoNorm(alpha-1,gamma)"
statpar=c(alpha-1,gamma)
rndpar1=c(nstart,1,statpar) #parámetros del generador de números aleatorios estacionario
stattimes=eval(statdist,rndpar1) #feval??? #tiempos de servicio estacionarios
arrtimes=rep(0,length(stattimes)) #size???
}
else {stattimes=NULL;arrtimes=NULL} #[]???
npoints=rpois(1,lambda*tmax) ##generar los arribos de los nuevos clientes según distrib Poisson
if(npoints>0) {arrtimes=rbind(arrtimes,sort(runif(npoints,1)*tmax))} #se condiciona que el número de puntos sea N, los puntos tienen distrib uniforme
ntotal=length(arrtimes)
 
rndpar2=c(ntotal-nstart,1,servpar) #parámetros para el generador de números aleatorios
servtimes=rbind(stattimes,feval(servdist,rndpar2)) #feval??? #genera los tiempos de servicio apra los nuevos clientes
deptimes=arrtime+servtimes #tiempos de salida
 
col2=c(rep(1,ntotal),rep(-1,ntotal)); col1=c(arrtimes,deptimes);arrate=data.frame(col1,col2)
arrate=with(arrate,arrate[order(arrate[,1]),]) #saltos ordenados
jumptimes=arrate[,1]
size=cumsum(arrate[,2])
 
ci=which(jumptimes>0 && jumptimes<=tmax) #cortar el proceso de conteo para que no exceda el tmax y borrar los tiempos 0 en el comienzo de los arribos "negativos"
if(ci>0) #~isempty(ci)??
{
fc=max(ci[1]-1)
lc=ci[length(ci)]
jumptimes=jumptimes[fc:lc]
size=size[fc:lc]
jumptimes=c(jumptimes,tmax)
size=c(size,size[length(size)])
}
plot(jumpstimes,type="o",pch=20)
plot(jumpstimes,size,type="o",pch=20)
out=list(jumpstimes=jumpstimes,size=size)
out #jumptimes: tiempos de cambios de estado en el sistema, size: número de clientes en el sistema
}
 
####Generación de procesos de nacimiento y muerte
 
#1)Intensidades de nacimiento y muerte racionales (lambda y mu fijas)
 
#En matlab: genera las intensidades de salto y tiempo de espera hasta el próximo salto
#lambda_i=lambda/(1+i); if i==0 mu_i=0; else mu_i=mu*i; end
#time=-log(rand)./(lambda_i+mu_i);
TBirthDeath<-function(lambda,mu,k) #en el nivel i la intensidad es lambda_i=lambda/(1+i) y la intensidad de muerte es mu_i=mu*i (siendo lambda y mu constantes fijadas)
{
time=rep(0,k)
Lambda=rep(0,k+1);Lambda[1]=lambda;
Mu=rep(0,k+1);Mu[1]=mu
for(i in 1:k)
{
lambda_i=lambda/(1+i); Lambda[i+1]=lambda_i
if(i==0) {mu_i=0} else {mu_i=mu*i}; Mu[i+1]=mu_i #mu*i??? gran valor!
time[i]=-log(runif(1))/(lambda_i+mu_i)
}
nf <- layout(matrix(c(1,2,3,3),2,2,byrow=TRUE), TRUE)
layout.show(nf)
plot(Lambda,main="Birth intensity",pch=20,type="o")
plot(Mu,main="Death intensity",pch=20,type="o")
plot(time,main="Time",pch=20,type="o")
out=list(Lambda=Lambda,Mu=Mu,time=time)
out
}
 
#BirthDeath: genera una trayectoria de procesos de nacimientos y muertes
BirthDeath=function(npoints=10,lambda=.5,mu=.2,grafico=FALSE) #npoints:largo de la trayectoria, flambda: la función-inline para calcular la intensidad de nacimientos en un punto dado (opcional), fmu: la función inline para calcular la intensidad de mortandad en un punto dado (opcional)
{
i=0 #comienza en el nivel i
tjump=rep(0,npoints); tjump[1]=0 #comienza en el tiempo 0
state=rep(0,npoints);state[1]=i #en el tiempo 0: nivel i
if(grafico==TRUE){Lambda=rep(0,npoints+1);Lambda[1]=lambda;Mu=rep(0,npoints+1);Mu[1]=mu}
for(k in 1:npoints)
{
lambda_i=lambda/(4+2*i);mu_i=4*mu*i #calcular las intensidades; Aquí: lambda=5 y mu=1???
time=-log(runif(1))/(lambda_i+mu_i) #tiempos entre saltos, con distribución exponencial de parámetro lambda_i+mu_i
if(runif(1)<=lambda_i/(lambda_i+mu_i)){i=i+1} else {i=i-1} #nacen o mueren, respectivamente
state[k]=i
tjump[k]=time
if(grafico==TRUE){Lambda[k+1]=lambda_i;Mu[k+1]=mu_i}
}
if(grafico==TRUE)
{
nf <- layout(matrix(c(1,2,3,3),2,2,byrow=TRUE), TRUE)
layout.show(nf)
plot(Lambda,main="Birth intensity",pch=20,type="o")
plot(Mu,main="Death intensity",pch=20,type="o")
plot(tjump,state,pch=20,type="p")
}
else {plot(tjump,state,pch=20,type="p")}
out=list(tjump=tjump,state=state)
out #tjump: tiempos de salto, state: estados de la cadena de ;Arkov
}
 
#2)Modelo de Moran (espacio de estados finitos y fronteras absorbentes)
#En matlab:
# x = 1 : N + 1;
# lambda(x) = (x-1).*(1-(x-1)./N);
# mu(x)=(x-1).*(1-(x-1)./N);
# q(x) = lambda(x) + mu(x); #exponential rate for time between jumps
TMoran=function(N) #está ok???
{
lambda=rep(0,N+1);mu=rep(0,N+1);q=rep(0,N+1)
for(x in 1:(N+1))
{
lambda[x]=(x-1)*(1-(x-1)/N)
mu[x]=(x-1)*(1-(x-1)/N)
q[x]=lambda[x]+mu[x] #tasa exponencial para los saltos entre tiempos
}
nf <- layout(matrix(c(1,2,3,3),2,2,byrow=TRUE), TRUE)
layout.show(nf)
plot(lambda,main="Birth intensity",pch=20,type="o")
plot(mu,main="Death intensity",pch=20,type="o")
plot(q,pch=20,type="o")
q
}
 
#genera la trayectoria de un proceso de tipo Moran, Ej: tiene como salida el número de genes de un alelo tipo A1 en una población de individuos haploides que pueden existir como A1 o A2
Moran=function(popsize=30,initsize=10) #el tamaño poblacional es "popsize" y el número inicial de individuos con alelo tipo A1 es "initsize"
{
A1nmb=rep(0,popsize); A1nmb[1]=initsize
x=initsize
i=1
while(x>1 && x<popsize+1)
{
lambda = (x-1)*(1-(x-1)/popsize)
mu =(x-1)*(1-(x-1)/popsize)
if(lambda/(lambda+mu)>runif(1)){x=x+1;A1nmb[i]=x} else {x=x-1;A1nmb[i]=x}
i=i+1
}
nmbsteps=length(A1nmb)
y=A1nmb[1:(nmbsteps-1)]
rate=((y-1)*(1-(y-1)/popsize))+((y-1)*(1-(y-1)/popsize))
jumptimes=cumsum(-log(runif(nmbsteps-1))/rate)
jumptimes=c(0,jumptimes)
plot(jumptimes,A1nmb,pch=20,type="o")
out=list(A1nmb=A1nmb,jumptimes=jumptimes)
out #A1nmb: los números sucesivos de alelos tipo A1 en la población
}
 
 
####Generación de procesos de ramificación
#simulación del proceso de Galton-Watson
#en tiempo discreto,
#tiempos de vida de longitud unitaria y,
#cada muerte de un individuo da lugar al nacimiento de un número aleatorio de niños
 
#Función offspring(k), que otorga un vector de antecedentes de una población de tamaño k; Proceso de Galton-Watson
#En Matlab:
# p = [1=2 0 1=2];z = [cumsum(p)];n = length(p); #nmb of possible offspring
# offmu = dot(0 : n - 1, p); u1 = sort(rand(1, k)); #offspring mean
# for j = 1 : n
# u(j) = length(find(u1< z(j)));
# end
# u = diff([0 u]); nu = u *(0 : n - 1)´;
 
offspring=function(k=10,p=c(1/2,0,1/2))
{
z=cumsum(p)
n=length(p) #nmb de posibles progenies
offmu=(0:(n-1))%*%p #media de la progenie
u1=sort(runif(k))
u=rep(0,n)
for(j in 1:n)
{
u[j]=length(which(u1<z[j]))
}
#u=c(0,u) ##"diff([0 u])" derivadas????
offspring=u*(0:(n-1)) #revisar?
sum(offspring) #revisar?
}
 
#genera una trayectorias de un proceso de ramificación de Galton-Watson
#con distribución de progenie p=(p0,p1,...,pn) ej: p=[0.3 0.4 0.2 0.1]; % Examples of offspring probabilitites
#comienza con un tamaño poblacional "initsize" en el tiempo 0
GaltonWatson=function(nmbgen=10,initsize=40,p=c(1/2,0,1/2)) #nmbgen: número de generaciones, initsize: tamaño poblacional inicial,p: vector de probabilidades de progenie (p0,...,pn)
{
popsize=rep(0,nmbgen);popsize[1]=initsize
k=1
while(k<=nmbgen)
{
popsize[k+1]=sum(offspring(popsize[k],p)) #salen 3 nº??? ##ver terminar la func offspring
k=k+1
}
x=0:nmbgen
plot(x,popsize,pch=20,type="o")
popsize #popsize: los tamaños poblacionales sucesivos
}