Ehrenfest Diffusion Process Simulation Code in R
Post date: Feb 20, 2011 5:23:24 AM
## code for question 4.1.a (Ehrenfest Diffusion Process), STAT 515, Haoying Wang
rm(list = ls(all = TRUE)) # clean workspace
ssize = 10000 # sample size
m=25 # total number of particles
x1 = rep(NA, ssize) # initial number of particles in box A
x2 = rep(NA, ssize) # initial number of particles in box B
x3 = rep(NA,ssize) # check if the total number of particles always equals to m
x <- 1:25
x1[1]=5
x2[1]=m-x1[1]
x3[1]=m
for (i in 2:ssize)
{
# randomly choose a number in [1,25]
k=sample(x,1)
if (k<=x1[i-1]) # k is euqal or less than x1, choose a particle from A, transfer to B
{x1[i]=x1[i-1]-1
x2[i]=x2[i-1]+1}
if (k>x1[i-1]) # k is larger than x1, choose a particle from B, transfer to A
{x1[i]=x1[i-1]+1
x2[i]=x2[i-1]-1}
x3[i]=x1[i]+x2[i]
}
## calculate expected number of particles in A and B:
sum(x1)/ssize
sum(x2)/ssize
sum(x3)/ssize
hist(x1) # draw a histogram of x1
## end of code