plot.traj <- function(traj,term,ntraj,s0,s1,...) {
# Description: plot sample trajectories of the posterior moments
# Input:
#   - traj: sample trajectories as in 'traj' argument returned by 'read.forwsim'
#   - term: terminal decisions and posterior moments for sample trajectories as returned by 'read.forwsim'
#   - ntraj: number of trajectories to be plotted
#   - b: if specified, boundaries will be plotted as well. s0 and s1 need to be specified as well
#   - s0: value of the posterior log-SD at which both boundaries cross. The posterior mean at this point is given by b[1]
#   - s1: value of the posterior log-SD corresponding to the posterior means b[2],b[3]
#   - ...: other arguments to 'plot'. Note that arguments 'type','xlim' and 'ylim' are set automatically so specifying
#          them will cause an error
# Ouput: nothing.
# Side effects: plot sample trajectories of the posterior moments saved in the forward simulation

logs <- c(s0,max(traj$logs))
m <- c(term$b0[1],term$b0[1]+(term$b1[1]-term$b0[1])*max(traj$logs-s0)/(s1-s0))
plot(logs,m,type='l',xlim=c(min(s0,min(traj$logs)),max(traj$logs)),ylim=range(traj$m),...)
par(new=TRUE)
m <- c(term$b0[1],term$b0[1]+(term$b2[1]-term$b0[1])*max(traj$logs-s0)/(s1-s0))
plot(logs,m,type='l',xlim=c(min(s0,min(traj$logs)),max(traj$logs)),ylim=range(traj$m),...)
for (i in 0:(ntraj-1)) {
  lines(traj$logs[traj$trt==i],traj$m[traj$trt==i],lty=(i+1))
  if (term$term[term$sim==0 & term$trt==i] == 1) {
    points(term$logs[term$sim==0 & term$trt==i],term$m[term$sim==0 & term$trt==i],pch=15,cex=1.5)
    text(term$logs[term$sim==0 & term$trt==i],term$m[term$sim==0 & term$trt==i],'Phase III',pos=2)
  } else {
    points(term$logs[term$sim==0 & term$trt==i],term$m[term$sim==0 & term$trt==i],pch=16,cex=1.5)
  }
}

}
