Shading are under the Normal Distribution

A demonstration

Introduction

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 or scalar quantities
• mean the mean of your univariate normal distribution
• sd the mean of your -univariate normal distribution
• lower.tail logical; if TRUE (default), probabilities are \(P(X ≤ x)\) otherwise, \(P(X > x)\).
• poly.col colour used for shaded region (default = “skyblue”)
• add2plot logical; if FALSE (default) this function will create a new plot, otherwise it will add to an existing plot
• add.prob logical; if TRUE (default), the probabilities associated with the shaded region will appear in text
• add.q logical; if TRUE (default), the quantile associated with the shaded probability will appear on the x-axis
• add.legend logical; if TRUE, a legend with the mean and standard deviation will be plotted in a legend
• ... arguments to be passed to the plot function

References
Adapted from code from here. Used in my lab here

To see the function, click the “CODE” button on the right-hand-side our source it here.

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", line=1)
        axis(1, at=q[1], labels = "", col="skyblue", col.axis="skyblue")
        axis(1, at=q[1], labels = "", col="skyblue", col.axis="skyblue", line=0.5)
        
        axis(1, at=q[2], labels = q[2], col="skyblue", col.axis="skyblue", line=1)
        axis(1, at=q[2], labels = "", col="skyblue", col.axis="skyblue")
        axis(1, at=q[2], labels = "", col="skyblue", col.axis="skyblue", line=0.5)
            }
    
    
    
  } 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)
      }
      
      if (add.q) {
        axis(1, at=q, labels = "", col="skyblue", col.axis="skyblue")
        axis(1, at=q, labels = q, col="skyblue", col.axis="skyblue", line=1)
        axis(1, at=q, labels = "", col="skyblue", col.axis="skyblue", line=0.5)
      }
      
      # 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)
      }
      if (add.q) {
        axis(1, at=q, labels = "", col="skyblue", col.axis="skyblue")
        axis(1, at=q, labels = q, col="skyblue", col.axis="skyblue", line=1)
        axis(1, at=q, labels = "", col="skyblue", col.axis="skyblue", line=0.5)
      } 
    
      
    }  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)
  }
}

Examples

Lower tail probabilities

\(P(Z \leq z)\)

Shade in the region \(P(Z \leq z)\) where \(Z\sim N(0,1)\) is the standard normal \(N(0,1)\)

myz <- -0.2533
snorm(myz)

The plot above is a visualization of P(Z < -0.2533 ) = 0.4 for Z~N(0, 1).

\(P(X \leq x)\)

Shades in the region \(P(X \leq x)\) where \(X\sim N(\mu = \texttt{mean},\sigma = \texttt{mean})\)

x <- 7.1 ; mymean = 7 ; mysd = 0.1
snorm(x, mean = mymean , sd = mysd,  add.legend = TRUE )

The plot above is a visualization of P(X < 7.1 ) = 0.841 for X~N( 7 , 0.1 )

Upper tail

\(P(Z > z)\)

Shade in the region \(P(Z > z)\) where \(Z\sim N(0,1)\) is the standard normal \(N(0,1)\)

snorm(myz,  lower.tail=FALSE)

The plot above is a visualization of P(Z > -0.2533 ) = 0.4 for Z~N(0, 1).

\(P(X > x)\)

Shades in the region \(P(X > x)\) where \(X\sim N(\mu = \texttt{mean},\sigma = \texttt{mean})\)

snorm(x, mean = mymean , sd = mysd,  add.legend = TRUE , lower.tail = FALSE)

The plot above is a visualization of P(X > 7.1 ) = 0.159 for X~N( 7 , 0.1 )

Inner range of values

\(P(a < Z < b)\)

Shade in the region \(P(a < Z < b)\) where \(Z\sim N(0,1)\) is the standard normal \(N(0,1)\)

a <- -1.1
b <- 2.0
q <- c(a,b)
snorm(q,  lower.tail=FALSE)

The plot above is a visualization of P( -1.1 < Z < 2 ) = -0.842 for Z~N(0, 1).

\(P(a < X < b)\)

Shade in the region \(P(a < X < b)\) where \(X\sim N(\mu = \texttt{mean},\sigma = \texttt{mean})\)

a <- 6.9
b <- 7.1
q <- c(a,b)
snorm(q,  mymean, mysd, lower.tail=FALSE)

The plot above is a visualization of P( 6.9 < X < 7.1 ) = 0 for X~N( 7 , 0.1 )