Skip to content

Instantly share code, notes, and snippets.

@CnrLwlss
Last active December 11, 2015 18:38
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save CnrLwlss/4643064 to your computer and use it in GitHub Desktop.
Save CnrLwlss/4643064 to your computer and use it in GitHub Desktop.
Demonstrating a function for simulating constrained random walks. Like a discrete Brownian bridge. http://cnr.lwlss.net/ConstrainedRandomWalk
# Simulating n steps of a biased random walk
# starting from x0 and with decreasing steps
# occurring with probability theta
randomWalk=function(n=100,x0=0,theta=0.5){
x=x0
res=array(x0,dim=n+1)
unifs=runif(n)
for(i in 0:(n-1)){
if (unifs[i+1]<=theta){x=x-1}else{x=x+1}
res[i+2]=x
}
return(res)
}
# Simulating n steps of a constrained random walk
# starting from x0 with a target after n steps of xtarg
constrainedWalk=function(n=100,x0=0,xtarg=0){
# Note that unless (xtarg-x0) and n are either both odd
# or both even, target cannot be reached exactly in n steps
x=x0
res=array(x0,dim=n+1)
unifs=runif(n)
for(i in 0:(n-1)){
theta=(1-(xtarg-x)/(n-i))/2
if (unifs[i+1]<=theta){x=x-1}else{x=x+1}
res[i+2]=x
}
return(res)
}
# Some example simulations and plots
n=500
x0=0
xtarg=50
sampno=50
pdf("RandomWalkFigures.pdf",width=16,height=8)
op<-par(oma=c(1,0,0,0),mar=c(4.1,4.1,0,0))
# Unconstrained, unbiased random walk
ensemble=replicate(sampno,randomWalk(n,x0,0.5))
plot(NULL,xlim=c(0,n),ylim=range(ensemble[,1:sampno]),xlab="Step number",ylab="x",cex.axis=1,cex.lab=2)
for(i in 1:sampno) points(0:n,ensemble[,i],type="s",col=grey(0.8*i/sampno))
# Biased random walk
theta=0.4
ensemble_up=replicate(sampno,randomWalk(n,x0,theta))
ensemble_down=replicate(sampno,randomWalk(n,x0,1-theta))
plot(NULL,xlim=c(0,n),ylim=range(rbind(ensemble_up,ensemble_down)),xlab="Step number",ylab="x",cex.axis=1,cex.lab=2)
for(i in 1:(sampno/2)) points(0:n,ensemble_up[,i],type="s",col="black")
for(i in 1:(sampno/2)) points(0:n,ensemble_down[,i],type="s",col="red")
# Constrained random walk
ensemble=replicate(sampno,constrainedWalk(n,x0,xtarg))
plot(NULL,xlim=c(0,n),ylim=range(ensemble[,1:sampno]),xlab="Step number",ylab="x",cex.axis=1,cex.lab=2)
for(i in 1:sampno) points(0:n,ensemble[,i],type="s",col=grey(0.8*i/sampno))
par(op)
dev.off()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment