Blog_simulación_procesos
#####################Random walks###############
####Simple random walk
# "Start at 0, jump one step up with probability p and down with probability 1-p"
#in Matlab: p=0.5;
#y=[0 cumsum(2.*(rand(1,n-1)<=p)-1)]; % n steps
#plot([0:n-1],y);
simple_RW<-function(p,n) #para el caso de movimiento arriba-abajo
{
y=c(0,cumsum(2*runif(n-1)<=p)-1)
plot(y,type="l")
}
####Random step size random walk
#Take any zero mean distribution for step size, e.g. uniform.
#in Matlab: x=rand(1,n)-1/2;
#y=[0 cumsum(x)-1)];
#plot([0:n-1],y);
RandomStep_RW<-function(n)
{
y=c(0,cumsum(runif(n)-1/2)-1)
plot(y,type="l")
}
####Brownian motion
#This is the continuous analog of symmetric random walk, each increment W(s+t)-W(s) is Gaussian with distribution N(0, t) and increments over disjoint intervals are independent. It is typically simulated as an approximating random walk in discrete time.
#in Matlab: n=1000;
#dt=1;
#y=[0 cumsum(dt^0.5.*randn(1,n))]; % standard Brownian motion
#plot(0:n,y);
BrownianMotion<-function(n=1000,dt=1)
{
y=c(0,dt^.5*rnorm(n))
plot(cumsum(y), type="l")
}
BrownianMotion2<-function(x = 0, t0 = 0, T = 1, N = 100)
{
if (T <= t0)
stop("wrong times")
dt <- (T - t0)/N
t <- seq(t0, T, length = N + 1) #t[i] = t0 + (T-t0)*i/N, i in 0:N. t0=0 para el MB geométrico
X <- ts(cumsum(c(x, rnorm(N) * sqrt(dt))), start = t0, deltat = dt)
return(invisible(X))
}
###########Poisson processes##################
#Generate random events so that: (i) events occur independent of each other, (ii) two or more events can not occur at the same time point, (iii) the events occur with constant intensity. Number of events N(t) that occur from the time zero up to time t is Poisson distributed with expected value lambda*t. The counting process N(t) is a Poisson process. The successive times between arrivals is Exp(lambda) distributed.
####Fixed number of steps
#in Matlab: % n simulated interarrival times from Exp(lambda)
#interarr=[0 -log(rand(1, n))./lambda];
#stairs(cumsum(interarr), 0:n);
PoissonProcess_FNS<-function(n,lambda)
{
y=c(0,cumsum(-log(runif(n))/lambda))
plot(y,type="l")
}
####Fixed interval of time, one process
#Fix a time interval [0 tmax]. The total number of events in the interval is Poisson distributed with expected value lambda*tmax. Conditioned on the number of events that occured, the events are uniformly distributed in the interval.
#in Matlab: % the number of points is Poisson-distributed
#npoints = poissrnd(lambda*tmax); % conditioned that number of points is N,% the points are uniformly distributed
#if (npoints>0)
# arrt = [0; sort(rand(npoints, 1)*tmax);
#else
# arrt = 0;
#end
#stairs(arrt, 0:npoints); % plot the counting process
PoissonProcess_FIT_1<-function(tmax,lambda)
{
npoints=rpois(lambda=lambda*tmax,1)
if(npoints>0){arrt=c(0,sort(runif(npoints)*tmax))} else {arrt=0}
plot(arrt)
}
####Fixed interval of time, N processes
#It is complicated to vectorize the previous algorithm. A far from optimal possibility is to use the same algorithm as for renewal processes with exponentially distributed interrenewal times.
#in Matlab: % add zero to the arrival times for nicer plots
#tarr = zeros(1, nproc); % sums of exponentialy distributed interarrival times, % as matrix columns
#i = 1;
#while (min(tarr(i,:))<=tmax)
# tarr = [tarr; tarr(i, :)-log(rand(1, nproc))/lambda];
# i = i+1;
#end
#stairs(tarr, 0:size(tarr, 1)-1); % plot the counting processes
PoissonProcess_FIT_N<-function(nproc,tmax,lambda)
{
tarr=NULL
tarr=c(0,nproc)
i=1
while(min(tarr[i])<=tmax) {tarr=c(tarr,tarr[i]-log(runif(nproc))/lambda); i=i+1}
plot(tarr,type="l")
}
#######Higher-dimensional versions Llamar a: library(rgl)
### 2-dimensional random walk
#Plots the points (u(k),v(k)), k=1:n, in the (u,v)-plane, where (u(k)) and (v(k)) are one-dimensional random walks. The sample code plots four trajectories of the same walk, four different colors.
#in Matlab:
#n=100000;
#colorstr=['b' 'r' 'g' 'y'];
#for k=1:4
# z=2.*(rand(2,n)<0.5)-1;
# x=[zeros(1,2); cumsum(z')];
# col=colorstr(k);
# plot(x(:,1),x(:,2),col);
# hold on
#end
#grid
RW_2D<-function(n)
{
z1=2*(runif(n))-1 #xi=sample(c(1,-1),1)
z2=2*(runif(n))-1 #yi=sample(c(1,-1),1)
x=rbind(rep(0,2),cbind(z1,z2))
plot(x,type="l")
}
###3-dimensional random walk. Same as above, in 3-dimensional space.
#in Matlab:
#p=0.5;
#n=10000;
#colorstr=['b' 'r' 'g' 'y'];
#for k=1:4
# z=2.*(rand(3,n)<=p)-1;
# x=[zeros(1,3); cumsum(z')];
# col=colorstr(k);
# plot3(x(:,1),x(:,2),x(:,3),col);
# hold on
#end
#grid
RW_3D<-function(n) #solo haremos una trayectoria
{
X=NULL
for(i in 1:3)
{
z=cumsum(2*runif(n)-1)
X=cbind(X,z)
}
plot3d(X[,1],X[,2],X[,3],pch=20,col="red",size=4,xlab="x",ylab="y",zlab="z")
}
###3D Brownian motion
#in Matlab:
#npoints = 5000;
#dt = 1;
#bm = cumsum([zeros(1, 3); dt^0.5*randn(npoints-1, 3)]);
#figure(1);
#plot3(bm(:, 1), bm(:, 2), bm(:, 3), 'k');
#pcol = (bm-repmat(min(bm), npoints, 1))./ ...
# repmat(max(bm)-min(bm), npoints, 1);
#hold on;
#scatter3(bm(:, 1), bm(:, 2), bm(:, 3), ...
# 10, pcol, 'filled');
#grid on;
#hold off;
BrownianMotion_2D<-function(npoints=500,dt=1,D=2)
{
BM=NULL
for(i in 1:D)
{
bm=c(0,rnorm(npoints-1,sd=dt^0.5))
BM=cbind(BM,cumsum(bm))
}
plot(BM[,1],BM[,2],type="l",xlab="x",ylab="y",main="Brownian Motion 2D")
if(D==3){ plot3d(BM[,1],BM[,2],BM[,3],pch=20,col="red",size=4,xlab="x",ylab="y",zlab="z",main="Brownian Motion 3D")}
BM
}
BrownianMotion_otro2D<-function(n=1000,fun=rnorm,dt=1) {
x<-cumsum(fun(n,sd=dt^0.5))
y<-cumsum(fun(n,sd=dt^0.5))
#FALTA
plot(x,y,type="l",main="Brownian Motion 2D")
}
###Poisson points in 2D and 3D space
#This is the generic model for placing points in space randomly and independently. In any fixed set of space there will be a Poisson distributed number of points with intensity proportional to the volume of that set. The number of points in any two disjoint sets are independent.
#in Matlab: % Poisson points in the unit cube, % intensity lambda
#lambda=100;
#nmb=poissrnd(lambda)
#x=rand(1,nmb);
#y=rand(1,nmb);
#z=rand(1,nmb);
#grid
#scatter3(x,y,z,5,5.*rand(1,nmb));
PoissonPoints_3D<-function(lambda,tmax)
{
nmb=rpois(lambda=lambda*tmax,n=1)
x=runif(nmb)
y=runif(nmb)
z=runif(nmb)
plot3d(x,y,z,pch=20,col="red",size=4,xlab="x",ylab="y",zlab="z")
}
PoissonPoints_2D<-function(lambda=.5,tmax=100)
{
nmb=rpois(lambda=lambda*tmax,n=1) #nº de puntos
x=runif(nmb)
y=runif(nmb)
plot(x,y,pch=20,xlab="x",ylab="y",main="Poisson Process")
cbind(x,y)
}
###Proceso de Poisson no-homogéneo
PoissonProcess_NH<-function(npts=100,tc=-50) #condicionado a npts
{
spatpts <- matrix(0,npts,2)
inpts <- 1
while(inpts <= npts)
{
xpt <- runif(1)
ypt <- runif(1)
lambda<-exp(tc*(xpt*ypt))
if(lambda > runif(1))
{
spatpts[inpts,1:2] <- cbind(xpt,ypt)
inpts <- inpts+1
}
}
par(pty="s")
plot(spatpts[,1],spatpts[,2],pch=20,main="Poisson Process NH")
spatpts
}
######The M/M/1 model
#This is a single server buffer model in continuous time. Arrivals to the system are determined by a Poisson process with intensity lambda. The server needs an exponentially distributed random time for each customer, mean service time 1/mu. The resulting system size N(t), t>=0, is a Markov process in continuous time which evolves as follows. Start from N(0)=n_0. Wait an exponential time with intensity lambda+mu (intensity lambda if n_0=0), then jump upwards with probability lambda/(lambda+mu) and downwards with probability mu/(lambda+mu). Repeat.
#The change of dynamics at N=0 is handled with a short loop in the beginning:
#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(i,mu,lambda,k)
{
x=rep(0,k)
for(j in 1:k)
{
if(i==0){mutemp=0} else {mutemp=mu}
if(runif(1)<=lambda/(lambda+mutemp)) {i=i+1} else {i=i-1}
x[j]=i
}
plot(x,type="l",main="M/M/1")
x
}
####Birth-and-death processes
#General birth and death intensities
#As an example we take a model where on level i the birth intensity is lambda_i=lambda/(1+i) and the death intensity mu_i=mu*i, lambda and mu fixed constants. A loop is required so that jump intensities and the waiting time till next jump can be updated, i.e.
#in Matlab:
#lambda_i=lambda/(1+i); if i==0 mu_i=0; else mu_i=mu*i; end
#time=-log(rand)./(lambda_i+mu_i);
Birth_Death<-function(i,lambda,mu,k)
{
time=rep(0,k)
for(j in 1:k)
{
lambda_i=lambda/(1+i)
if(i==0) {mu_i=0} else {mu_i=mu*i}
time[j]=-log(runif(1))/(lambda_i+mu_i)
}
plot(time,type="l",main="Birth Death")
time
}
Comentarios
Publicar un comentario