# -------------------------------------------------------------------------------------------------------------------------  #
  #  -----------------------------------------------------------------------------------------------------------------------   #
  #                                                                                                                            #
  #           File Name              :  beta plot new.r                                                                        #
  #           Last Updated Funclist  :  15 Feb 2014,  1:13 AM (Saturday)                                                       #
  #                                                                                                                            #
  #           Author Name            :  Rick Saporta                                                                           #
  #           Author Email           :  RickSaporta@gmail.com                                                                  #
  #           Author URL             :  www.github.com/rsaporta                                                                #
  #                                                                                                                            #
  #           Packages Called        :  NA                                                                                     #
  #           Packages Used via NS   :  NA                                                                                     #
  #                                                                                                                            #
  #  -----------------------------------------------------------------------------------------------------------------------   #
  #                                                                                                                            #
  #   BetaPlot           ( a.seq, b.seq, x.min=0, x.max=1, pause=0.35, plot.this.many=length(a.seq), ci=0.9 )                  #
  #                                                                                                                            #
  #                                                                                                                            #
  #                                                     <END FUNCS>                                                            #
  #  -----------------------------------------------------------------------------------------------------------------------   #
  # -------------------------------------------------------------------------------------------------------------------------  #


BetaPlot <- function(a.seq, b.seq, x.min=0, x.max=1, pause=0.35, plot.this.many=length(a.seq), ci=.9)  {

  # color paramters
  color.first.plot  <- rgb(1, .10, .30, 1)
  color.last.plot   <- rgb(0, .80, 0,  1)
  color.scale.start <- rgb(.8, .1, .7, .05)
  color.scale.end   <- rgb(0, .8, .1, 0.9)
  color.mde         <- rgb(r=0, g=0.5, b=0.8, a=0.4)
  color.ci          <- rgb(r=0, g=0.5, b=0.9, a=0.8)


  x <- seq(from=x.min, to=x.max, length.out=250)

  # a & b should be the same length
  if (length(a.seq) != length(b.seq))
    stop("`a.seq` and `b.seq` must be the same length.")

  # caluclate Beta, then set names
  output <- mapply(Beta, a=a.seq, b=b.seq, MoreArgs=list(x=x), SIMPLIFY=FALSE)
  output <- setNames(output, paste0("a=", a.seq, ",  b=", b.seq))
  L <- length(output)

  # default to the length of the number of arguments, if 
  if (missing(plot.this.many) || !is.numeric(plot.this.many))
    plot.this.many <- L

  if (plot.this.many < L) {
    # if only 2, then select first and last
    if (plot.this.many==2)
      selected <- c(1, L)
    else 
      selected <- ceiling(c(1, seq(from=2, to=L, length.out=(plot.this.many-1))))
  
    # crop output and reset L. (note: output is already named)
    output <- output[selected]
    L <- length(output)
  }

# legend('topright', names(a)[-1] , 
#    lty=1, col=c('red', 'blue', 'green',' brown'), bty='n', cex=.75)


  # check for user error on the value of pause
  if (isTRUE(pause))
    pause <- 0.35
  if(!is.numeric(pause))
    pause <- 0
  if (pause < 0) {
    warning("pause should be a positive numeric value in seconds. Setting to 0.")
    pause <- 0
  }

  ##  -------------  PLOT THE MULTIPLE GRAPHS    -------------  ##
    # graph parameters
    y.max <- max(unlist(output), na.rm=TRUE)
    y.min <- min(unlist(output), na.rm=TRUE)
    x.min <- min(x, na.rm=TRUE)
    x.max <- max(x, na.rm=TRUE)
   
    colors <- make.rgba(L, start=color.scale.start, end=color.scale.end, color="log", final.alpha.gap=0.5)
    lwds   <- seq(0.2, 1.5, length.out=L)
    ltys    <- rep(4, L)

    # adjust first / last parameter
    colors[c(1, L)] <- c(color.first.plot, color.last.plot)
    lwds[c(1, L)]   <- 2
    ltys[c(1, L)]   <- 1

    # empty plot with correct x, y limits
    plot(x=x.min, y=y.min, ylim=c(y.min, y.max), xlim=c(x.min, x.max), type="n")
    
    # iterate
    for (i in 1:L) {
      Sys.sleep(pause)
      lines(x=x, y=output[[i]], lwd=lwds[[i]], lty=ltys[[i]], col=colors[[i]])
    }
  ##  -------------  PLOT THE MULTIPLE GRAPHS    -------------  ##


browser();browser();browser();

  ## find some stats
  out <- output[[L]]
  
  # mode
    mde.val  <- max(out) 
    mde.indx <- which(out == mde.val)
    mde.x.val <- x[[mde.indx]] 
    color.mde <- color.mde
    abline(v=mde.x.val, lwd=0.5, lty=8, col=color.mde)

  # ci bounds
    ci.lims <- find.ci.single.mode(out, ci=ci, use="mode")
    x.intv  <- x[ci.lims]
    abline(v=x[ci.lims], lwd=1, lty=2, col=color.ci)

    # horizontal bar
    lines(x=x[ci.lims], y=out[ci.lims], lwd=1.5, lty=2, col=color.ci)

    names(x.intv) <- names(ci.lims)
    x.intv <-  fw3(x.intv, digs=2)
    md <- paste0("Posterior Mode is at ", round(mde.x.val, 3))
    ci.msg <- paste0("A ", fwp(ci, 0, sep=""), " credible interval\n"
                    ,"on the posterior is ["
                    , x.intv[["lower"]], ", ", x.intv[["upper"]], "]")

  legend("topleft", c(names(output), md))
  title(paste("Beta distribution", ci.msg, sep="\n"))
  return(output[[L]])
}
BetaPlot(seqs$a, seqs$b, plot.this.many=10, pause=0.1)