RW = function (n, M, p = 0.5) 
{
    pos <- n
    rw <- pos
    status <- 0
    if (pos == M) 
        status <- "win"
    if (pos == 0) 
        status <- "ruin"
    while (status == 0) {
        x <- rbinom(1, 1, p)
        if (x == 1) 
            pos <- pos + 1
        if (x == 0) 
            pos <- pos - 1
        rw <- c(rw, pos)
        if (pos == M) 
            status <- "win"
        if (pos == 0) 
            status <- "ruin"
    }
    l <- list(rw = rw, result = status)
    return(l)
}


simRW = function(num, n, M, p=0.5, PLOT=TRUE, steps=-1) {
  # num is the number of trials; 
  # steps is the number of steps to display on the plot
  # if "steps" is omitted, the plot will resize automatically to the length of the longest walk
  # if PLOT is TRUE, the plot is output, otherwise not
  
  q <- 1 - p
  
  cat(paste("Gambler's ruin, p=", p, ", n=", n, ", M=", M, "\n",sep=""))
  
  if (p==0.5) ruinprob <- n/M
  else 
  {
    ruinprob <- 1-   (  (1-(q/p)^n)    /    (1-(q/p)^M) )
  }
  cat(paste("Probability to be absorbed at 0 is ", round(ruinprob,5), "\n", sep=""))
  
  if (p==0.5) avtime <- n*(M-n)
  else
  {
    avtime <- M*(1-(q/p)^n)/((p-q)*(1-(q/p)^M))-n/(p-q)
  }
  cat(paste("Expected run length (until absorption at 0 or M) is ", round(avtime,2), "\n", sep=""))
  
  seed=sample(1:1000000,1); set.seed(seed)
  # we set a random seed and remember it since we will run the random walks twice, 
  # once to find the longest length and the second time to plot them all
  
  maxlength <- 0  
	count <- 0
	exp.length = rep(0, num)
	for(i in 1:num) 
    {
          a  = RW(n, M, p)
          exp.length[i] = (length(a$rw) - 1)
          if(length(a$rw)>maxlength) maxlength = length(a$rw)
          if(a$result == "ruin") count = count + 1
	   }
  cat(paste(count, " out of ", num, " trials absorbed at 0\n", sep=""))
	cat(paste("Proportion absorbed at 0: ", round(count/num,5), "\n", sep = ""))
	cat(paste("Average observed run length: ", mean(exp.length), "\n", sep = ""))

  set.seed(seed)

  if(PLOT)
    {  
    if(steps>0) displaylength=steps else displaylength=maxlength
  
    for(i in 1:num) {
      a  = RW(n, M, p)
      if(i == 1) {
        plot(1:length(a$rw) - 1, a$rw, xlim = c(0, displaylength), ylim = c(0, M), typ = "s",
             lwd = 3, pch=20, xlab = "Number of Steps", ylab = "Position", 
             main =  paste("Gambler's Ruin: p = ", round(p, 2), 
                            ", n = ",round(n, 0), 
                              ", M = ", round(M, 0), sep = "")) } 
      else {
        lines(  1:length(a$rw) - 1, a$rw, col = i, lwd = 3, typ="s", pch=20)
      }
    }
  }
}