How to plot a normal distribution by labeling specific parts of the x-axis?

An extremely inefficient and unusual, but beautiful solution, which works based on the ideas of Monte Carlo simulation, is this:

1. simulate many draws (or samples) from a given distribution (say the normal).
2. plot the density of these draws using rnorm. The rnorm function takes as arguments (A,B,C) and returns a vector of A samples from a normal distribution centered at B, with standard deviation C.

Thus to take a sample of size 50,000 from a standard normal (i.e, a normal with mean 0 and standard deviation 1), and plot its density, we do the following:

x = rnorm(50000,0,1)

plot(density(x))

As the number of draws goes to infinity this will converge in distribution to the normal. To illustrate this, see the image below which shows from left to right and top to bottom 5000,50000,500000, and 5 million samples.

The easiest (but not general) way is to restrict the limits of the x axis. The +/- 1:3 sigma will be labeled as such, and the mean will be labeled as 0 - indicating 0 deviations from the mean.

plot(x,y, type = "l", lwd = 2, xlim = c(-3.5,3.5))

Another option is to use more specific labels:

plot(x,y, type = "l", lwd = 2, axes = FALSE, xlab = "", ylab = "")
axis(1, at = -3:3, labels = c("-3s", "-2s", "-1s", "mean", "1s", "2s", "3s"))

Using the code in this answer, you could skip creating x and just use curve() on the dnorm function:

curve(dnorm, -3.5, 3.5, lwd=2, axes = FALSE, xlab = "", ylab = "")
axis(1, at = -3:3, labels = c("-3s", "-2s", "-1s", "mean", "1s", "2s", "3s"))

But this doesn't use the given code anymore.

If you like hard way of doing something without using R built in function or you want to do this outside R, you can use the following formula.

x<-seq(-4,4,length=200)
s = 1
mu = 0
y <- (1/(s * sqrt(2*pi))) * exp(-((x-mu)^2)/(2*s^2))
plot(x,y, type="l", lwd=2, col = "blue", xlim = c(-3.5,3.5))