snorm for shaded normal density curve

Adapted from code from here. Used in my lab here

Description
Plots the pdf of the normal distribution with mean equal to mean and standard deviation equal to sd and shades in the area under the curve as determined by q (either a scalar of a vector of length 2)

Usage

qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, ...)

• q vector of quantiles

snorm <- function(q, mean=0, sd=1, rsd=4, lower.tail=TRUE, poly.col="skyblue", 
                  add2plot = FALSE, add.prob=TRUE, add.q=TRUE, add.legend=FALSE, ...){
  
  if (length(q)==2) {
    q <- sort(q)
    xx=seq(mean-rsd*sd,mean+rsd*sd,length=200)
    yy=dnorm(xx,mean,sd)
    if (!add2plot) {
      plot(xx,yy,type="l", xlab="x", ylab=  expression(paste(f*"("* x *" ; " * mu *", "*sigma*")")), ...)
    } else {
      lines(xx, yy, type="l")
    }
    x=seq(q[1],q[2],length=100)
    y=dnorm(x,mean,sd)
    polygon(c(q[1],x,q[2]),c(0,y,0),col=poly.col)
    
    if (add.prob){
      textheight <- max(dnorm(x, mean, sd))/4
      prob <- round(pnorm(q[2], mean, sd) - pnorm(q[1], mean, sd), 3)
      text(mean(q),textheight, prob)
    }
    if (add.q){
      axis(1, at=q[1], labels = q[1], col="skyblue", col.axis="skyblue")
      axis(1, at=q[2], labels = q[2], col="skyblue", col.axis="skyblue")
    }
    
  } else {
    if (lower.tail) {
      # lower tail probabilities:
      xx=seq(mean-rsd*sd,mean+rsd*sd,length=200)
      yy=dnorm(xx,mean,sd)
      if (!add2plot) {
        plot(xx,yy,type="l", xlab="x", ylab=  expression(paste(f*"("* x *" ; " * mu *", "*sigma*")")) , ...)
      }  else {
        lines(xx, yy, type="l")
      }
      
      x=seq(xx[1],q,length=100)
      y=dnorm(x,mean, sd)
      polygon(c(xx[1],x,q),c(0,y,0),col=poly.col)
      if (add.prob){
        textheight <- dnorm(q, mean, sd)/4
        prob <- round(pnorm(q, mean, sd), 3)
        text(q-0.5*sd,textheight, prob)
      }
      
      # arrows(0.5,0.1,-0.2,0,length=.15)
      # text(0.5,0.12,"-0.2533")
      
    } else if (!lower.tail) {
      
      xx=seq(mean-rsd*sd,mean+rsd*sd,length=200)
      yy=dnorm(xx,mean,sd)
      if (!add2plot) {
        plot(xx,yy,type="l", xlab="x", ylab=  expression(paste(f*"("* x *" ; " * mu *", "*sigma*")")), ...)
      } else {
        lines(xx, yy, type="l")
      }
      
      x=seq(q,xx[length(xx)],length=100)
      y=dnorm(x,mean,sd)
      polygon(c(q,x,xx[length(xx)]),c(0,y,0),col=poly.col)
      if (add.prob){
        textheight <- dnorm(q, mean, sd)/4
        prob <- round(pnorm(q, mean, sd, lower.tail = FALSE), 3)
        text(q+0.5*sd,textheight, prob)
      }
      
    }  else {
      stop("please specify an appropriate value for lower.tail: TRUE, FALSE, or inner")
    }
  }
  

  if (add.legend)  {
    ltext <- bquote(mu *" = "* .(mean)*", "*sigma *" = "* .(sd))
    legend("topright", lwd= 2,
       legend = ltext)
  }
}

snorm(-0.2533)

snorm(0.5244, lower.tail=FALSE)