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