\begin{Schunk}
\begin{Sinput}
 ###Set mu and sigma at your pleasure:
 mu <- 10.03
 sigma <- 12.5786786
 myx <- seq( mu - 3.5*sigma,  mu+ 3.5*sigma, length.out=500)
 myDensity <- dnorm(myx,mean=mu,sd=sigma)
 # It is challenging to combine plotmath with values of mu and sigma in one expression.
 # Either bquote or substitute can be used.  First use bquote:
 
 myTitle1 <- bquote (paste("x ~ Normal(", mu== .(round(mu,2)), ',', sigma== .(round(sigma,2)),")") )
 ### Using substitute:
 ### myTitle1 <-  substitute( "x ~ Normal" ~~ group( "(", list(mu==mu1, sigma^2==sigma2)#, ")") ,  list(mu1=round(mu,2), sigma2=round(sigma^2,2)))
 
 ### Or combine the two:
 ### t1 <- substitute ( mu == a ,  list (a = mu))
 ### t2 <- substitute ( sigma == a, list(a = round(sigma,2)))
 ### myTitle1 <- bquote(paste("x ~ Normal(", .(t1),",", .(t2),")" ) )
 ### xpd needed to allow writing outside strict box of graph
 par(xpd=TRUE, ps=10)
 plot(myx, myDensity, type="l", xlab="x", ylab="Probability Density ", main=myTitle1, axes=FALSE)
 axis(2, pos= mu - 3.6*sigma)
 axis(1, pos=0)
 lines(c(myx[1],myx[length(myx)]),c(0,0)) ### closes off axes
 # bquote creates an expression that text plotters can use
 t1 <-  bquote( mu== .(mu))
 # Find a value of myx that is "very close to" mu
 centerX <- max(which (myx <= mu))
 # plot light vertical line under peak of density
 lines( c(mu, mu), c(0, myDensity[centerX]), lty= 14, lwd=.2)
 # label the mean in the bottom margin
 mtext(bquote( mu == .(mu)), 1, at=mu, line=-1)
 ### find position 20% "up" vertically, to use for arrow coordinate
 ss = 0.2 * max(myDensity)
 # Insert interval to represent width of one sigma
 arrows( x0=mu, y0= ss, x1=mu+sigma, y1=ss, code=3, angle=90, length=0.1)
 ### Write the value of sigma above that interval
 t2 <-  bquote( sigma== .(round(sigma,2)))
 text( mu+0.5*sigma, 1.15*ss, t2)
 ### Create a formula for the Normal
 normalFormula <- expression (f(x) == frac (1, sigma* sqrt(2*pi)) * e^{~~ - ~~ frac(1,2)~~ bgroup("(", frac(x-mu,sigma),")")^2})
 # Draw the Normal formula
 text ( mu + 0.5*sigma, max(myDensity)- 0.10 * max(myDensity),  normalFormula, pos=4)
 ### Theory says we should have 2.5% of the area to the left of: -1.96 * sigma.
 ### Find the X coordinate of that "critical value"
 criticalValue <- mu -1.96 * sigma
 criticalValue <- round(criticalValue, 2)
 ### Then grab all myx values that are "to the left" of that critical value.
 specialX <-  myx[myx <= criticalValue]
 ### mark the critical value in the graph
 text ( criticalValue, 0 ,
       label=criticalValue, pos=1)
 ### Take sequence parallel to values of myx inside critical region
 specialY <- myDensity[myx < criticalValue]
 #  Polygon makes a nice shaded illustration
 polygon(c(specialX[1], specialX, specialX[length(specialX )]), c(0, specialY, 0), density=c(-110),col="lightgray" )
 ### I want to insert annotation about area on left side.
 shadedArea <- round(pnorm(mu - 1.96 * sigma, mean=mu, sd=sigma),4)
 al1 <- bquote(Prob(x <= .(criticalValue)))
 al2 <- bquote(phantom(0) == F( .(criticalValue) ))
 al3 <- bquote(phantom(0) == .(round(shadedArea,3)))
 ### Hard to position this text "just right"
 ### Have tried many ideas, this may be least bad.
 ### Get center position in shaded area
 medX <- median(specialX)
 ### density at that center point:
 denAtMedX <- myDensity[indexMed <- max(which(specialX < medX))]
 text(medX, denAtMedX+0.0055, labels = al1)
 text(medX, denAtMedX+0.004, labels = al2)
 text(medX, denAtMedX+0.0025, labels = al3)
 ### point from text toward shaded area
 arrows( x0 = medX, y0 = myDensity[indexMed]+0.002,
        x1 = mu-2.5 *sigma, y1 = 0.40*myDensity[length(specialX)], length = 0.1)
 ss <- 0.1 * max(myDensity)
 ### Mark interval from mu to mu-1.96*sigma
 arrows( x0 = mu, y0 = ss, x1 = mu-1.96*sigma, y1 = ss, code = 3, angle = 90, length = 0.1)
 ### Put text above interval
 text( mu - 2.0*sigma,1.15*ss, bquote(paste(.(criticalValue), phantom(1)==mu-1.96 * sigma,sep="")), pos=4)
 criticalValue <- mu +1.96 * sigma
 criticalValue <- round(criticalValue, 2)
 ### Then grab all myx values that are "to the left" of that critical value.
 specialX <-  myx[myx >= criticalValue]
 ### mark the critical value in the graph
 text ( criticalValue, 0 , label= paste(criticalValue), pos=1)
 ### Take sequence parallel to values of myx inside critical region
 specialY <- myDensity[myx > criticalValue]
 #  Polygon makes a nice shaded illustration
 polygon(c(specialX[1], specialX, specialX[length(specialX )]), c(0, specialY, 0), density=c(-110), col = "lightgray" )
 shadedArea <- round(pnorm(mu + 1.96 * sigma, mean = mu, sd = sigma, lower.tail = F),4)
 ### Insert simpler comment on right side.
 
 al2 <- bquote(phantom(0) == 1 - F( .(criticalValue) ))
 al3 <- bquote(phantom(0) == .(round(shadedArea,3)))
 ### Hard to position this text "just right"
 ### Have tried many ideas, this may be least bad.
 ### Get center position in shaded area
 medX <- median(specialX)
 ### density at that center point:
 denAtMedX <- myDensity[indexMed <- max(which(specialX < medX))]
 text(medX, denAtMedX+0.004, labels = al2)
 text(medX, denAtMedX+0.0025, labels = al3)
 ### point from text toward shaded area
 arrows(x0 = medX, y0 = myDensity[indexMed]+0.002, x1 = mu+2.5 *sigma, y1 = 0.40*myDensity[length(specialX)] ,   length = 0.1)
\end{Sinput}
\end{Schunk}
