gnorm <- function(mu, sigma,prob=F,a=NA,b=NA,quantile=F,p=0.5)
{
  values <- seq(-1,1,.005) * 4 * sigma + mu
  probs <- dnorm(values, mu, sigma)
  plot(values, probs, axes = F, type = "n", xlab = "Possible Values", 
    ylab = "Probability Density", main = paste("Normal Distribution \n", 
    "mu =", mu, ", sigma =", sigma))
  axis(1, pos = 0)
  abline(0,0,col=1)
  lines(values, probs, col = 2)
  lo <- mu - 4 * sigma
  hi <- mu + 4 * sigma
  if(prob==T)
  {
    if(!is.na(a) && !is.na(b) && a > b){
      d <- a
      a <- b
      b <- d
    }
    if(is.na(a) || a <= lo)
      ulo <- lo
    else if(a <= hi)
      ulo <- a
    else
      ulo <- hi
    if(is.na(b) || b >= hi)
      uhi <- hi
    else if(b >= lo)
      uhi <- b
    else
      uhi <- lo
#    if(!is.na(a) && !is.na(b)){
#      if(a > lo && b < hi)
#      {
#	 u <- seq(a,b,length=601)
#      }
#      else if(a <= lo && b >= lo && b < hi)
#      {
#	 u <- seq(lo,b,length=601)
#      }
#      else if(a > lo && a <= hi && b >= hi)
#      {
#	 u <- seq(a,hi,length=601)
#      }
#      else if(a <= lo && b >= hi)
#      {
#	 u <- seq(lo,hi,length=601)
#      }
#    }
    u <- seq(ulo,uhi,length=601)
    lines(u,dnorm(u,mu,sigma),type="h",col=2)
    if(!is.na(a) && !is.na(b)){
      text(mu - 3.9 * sigma, 0.8 * dnorm(mu,mu,sigma),
        paste("P( ",a," < X < ",b," ) = ",
  	round(pnorm(b,mu,sigma)-pnorm(a,mu,sigma),digits=4),sep=""),
        adj=0,col=4,cex=1.2)
      text(mu - 3.9 * sigma, 0.5 * dnorm(mu,mu,sigma),
        paste("P( X < ",a," ) = ",
          round(pnorm(a,mu,sigma),digits=4),sep=""),adj=0,col=4,cex=1.2)
      text(mu + 3.9 * sigma, 0.5 * dnorm(mu,mu,sigma),
        paste("P( X > ",b," ) = ",
       	round(1-pnorm(b,mu,sigma),digits=4),sep=""),adj=1,col=4,cex=1.2)
    }
    else if(!is.na(a) && is.na(b)){
      text(mu - 3.9 * sigma, 0.5 * dnorm(mu,mu,sigma),
        paste("P( X < ",a," ) = ",
          round(pnorm(a,mu,sigma),digits=4),sep=""),adj=0,col=4,cex=1.2)
      text(mu + 3.9 * sigma, 0.5 * dnorm(mu,mu,sigma),
        paste("P( X > ",a," ) = ",
       	round(1-pnorm(a,mu,sigma),digits=4),sep=""),adj=1,col=4,cex=1.2)
    }
    else if(is.na(a) && !is.na(b)){
      text(mu - 3.9 * sigma, 0.5 * dnorm(mu,mu,sigma),
        paste("P( X < ",b," ) = ",
          round(pnorm(b,mu,sigma),digits=4),sep=""),adj=0,col=4,cex=1.2)
      text(mu + 3.9 * sigma, 0.5 * dnorm(mu,mu,sigma),
        paste("P( X > ",b," ) = ",
       	round(1-pnorm(b,mu,sigma),digits=4),sep=""),adj=1,col=4,cex=1.2)
    }
  }
  else if(quantile==T)
  {
    if( p <= 0 || p >= 1) p <- 0.5
    x <- qnorm(p,mu,sigma)
    if( x > lo && x < hi)
    {
      u <- seq(lo,x,length=601)
      lines(u,dnorm(u,mu,sigma),type="h",col=2)
      text(x, -0.09 * dnorm(mu,mu,sigma),
  	paste("z = ",round(qnorm(p),2),sep=""),adj=0.5,col=4)
    }
    else if(x >= hi)
    {
      u <- seq(lo,hi,length=601)
      lines(u,dnorm(u,mu,sigma),type="h",col=2)
      text(hi, -0.09 * dnorm(mu,mu,sigma),
  	paste("z = ",round(qnorm(p),2),sep=""),adj=0.5,col=4)
    }
    else if( x <= lo)
    {
      text(lo, -0.09 * dnorm(mu,mu,sigma),
  	paste("z = ",round(qnorm(p),2),sep=""),adj=0.5,col=4)
    }
    text(mu - 3.9 * sigma, 0.3 * dnorm(mu,mu,sigma),
      paste("P( X < ",signif(x,4)," ) = ",
  	round(p,digits=4),sep=""),adj=0,col=4)
  }
  return(invisible())
}

gbinom <- function(n, p, low=0, high=n,scale = F)
{
  sd <- sqrt(n * p * (1 - p))
  if(scale && (n > 10)) {
	low <- max(0, round(n * p - 4 * sd))
	high <- min(n, round(n * p + 4 * sd))
  }
  values <- low:high
  probs <- dbinom(values, n, p)
  plot(c(low,high), c(0,max(probs)), type = "n", xlab = "Possible Values", ylab = 
    "Probability", main = paste("Binomial Distribution \n", "n =", 
    n, ", p =", p))
  lines(values, probs, type = "h", col = 2)
  abline(h=0,col=3)
  return(invisible())
}