"poisson.process"<-
function(lambda = 0.1, N = 30, initial.size = 10)
{
	#
	# To set the random number seed, un-comment the following line 
	# and plug in any non-negative integer in [0,1023].
	# set.seed(123)
	#
	# Simulate N events from a Poisson process.
	# 
	S <- rep(0, N)
	W <- S
	X <- S
	X[1] <- initial.size
	W[1] <- S[1]
	#
	# Note that this is time = 0.
	#
	for(i in 2:N) {
		#
		# Simulate a sojourn (inter-event) time, and compute the waiting time.
		#
		S[i] <- rexp(1, rate = lambda)
		W[i] <- W[i - 1] + S[i]
		X[i] <- X[i - 1] + 1
	}
	# end for loop on i
	#
	# Compute the analytic mean function of the process on a fine grid of times.
	# 
	grid <- ((0:1000) * max(W))/1000
	M <- initial.size + (lambda * grid)
	#
	# Set up a plotting region, then add curve for the analytical mean and step function
	# for the Poisson process.
	#
	plot(c(W, grid), c(X, M), type = "n", xlab = "Time", ylab = 
		"Population Size")
	lines(grid, M, col = 2)  
	#col = 2 is red in R
	# col = 8 is red in Splus
	lines(W,X,type="s")
	return()
}