library(tidyverse)   # Only needed for the plot below. The probability functions are Core R.

# Example - P(X≤6) for a N(3,2^2) Distribution
########################################################################################

# Plot with the P(X≤6) area shaded
x_axis <- tibble(
  xlim = c(-4, 10)
)
ggplot(x_axis , aes(x = xlim)) +
  stat_function(fun = dnorm, args = list(mean = 3, sd = 2)) +
  stat_function(fun = dnorm, args = list(mean = 3, sd = 2),
                geom = "area", fill = "steelblue", alpha = .25, xlim=c(-4,6)) +
  labs(x = "X", y = "Density")

pnorm(6, mean=3, sd=2)    # P(X ≤ 6) for N(3,2^2)
pnorm(1.5, mean=0, sd=1)  # P(Z ≤ 1.5) for N(0,1) Standard Normal, using the z-score (6-3)/2 = 1.5

# Symmetry: 0 and 6 are both 3 units from the mean, so P(X ≥ 0) = P(X ≤ 6)
pnorm(0, mean=3, sd=2, lower.tail=FALSE)   # P(X ≥ 0)
1 - pnorm(0, mean=3, sd=2)                 # P(X ≥ 0) = 1 - P(X ≤ 0)


# Example - 68-95-99.7 Rule
########################################################################################

# Calculations for a N(0,1) standard normal distribution
1-2*pnorm(-1, mean=0, sd=1)  # ±1 sd on each side of mean
1-2*pnorm(-2, mean=0, sd=1)  # ±2 sd on each side of mean
1-2*pnorm(-3, mean=0, sd=1)  # ±3 sd on each side of mean

# Calculations for a N(10,5^2) normal distribution
pnorm(15, mean=10, sd=5) - pnorm(5, mean=10, sd=5)  # ±1 sd on each side of mean
pnorm(20, mean=10, sd=5) - pnorm(0, mean=10, sd=5)  # ±2 sd on each side of mean
pnorm(25, mean=10, sd=5) - pnorm(-5, mean=10, sd=5) # ±3 sd on each side of mean