--- title: "Introduction to dnorm, pnorm, qnorm, and rnorm for new biostatisticians" author: "Sean Kross" date: "October 1, 2015" output: html_document --- Today I was in Dan's office hours and someone asked, "what is the equivalent in R of the back of the stats textbook table of probabilities and their corresponding Z-scores?" ([This](http://home.ubalt.edu/ntsbarsh/Business-stat/StatistialTables.pdf) is an example of the kind of table the student was talking about.) This question indicated to me that although we've been asked to use some of the distribution functions in past homeworks, there may be some misunderstanding about how these functions work. Right now I'm going to focus on the functions for the normal distribution, but you can find a list of all distribution functions by typing `help(Distributions)` into your R console. --- ## `dnorm` As we all know the probability density for the normal distribution is: $$ f(x|\mu,\sigma) = \frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-\mu)^2}{2\sigma^2}} $$ The function `dnorm` returns the value of the probability density function for the normal distribution given parameters for $x$, $\mu$, and $\sigma$. Some examples of using `dnorm` are below: ```{r} # This is a comment. Anything I write after the octothorpe is not executed. # This is the same as computing the pdf of the normal with x = 0, mu = 0 and # sigma = 0. The dnorm function takes three main arguments, as do all of the # *norm functions in R. dnorm(0, mean = 0, sd = 1) # The line of code below does the same thing as the same as the line of code # above, since mean = 0 and sd = 0 are the default arguments for the dnorm # function. dnorm(0) # Another exmaple of dnorm where parameters have been changed. dnorm(2, mean = 5, sd = 3) ``` Although $x$ represents the independent variable of the pdf for the normal distribution, it's also useful to think of $x$ as a Z-score. Let me show you what I mean by graphing the pdf of the normal distribution with `dnorm`. ```{r} # First I'll make a vector of Z-scores z_scores <- seq(-3, 3, by = .1) # Let's print the vector z_scores # Let's make a vector of the values the function takes given those Z-scores. # Remember for dnorm the default value for mean is 0 and for sd is 1. dvalues <- dnorm(z_scores) # Let's examine those values dvalues # Now we'll plot these values plot(dvalues, # Plot where y = values and x = index of the value in the vector xaxt = "n", # Don't label the x-axis type = "l", # Make it a line plot main = "pdf of the Standard Normal", xlab= "Z-score") # These commands label the x-axis axis(1, at=which(dvalues == dnorm(0)), labels=c(0)) axis(1, at=which(dvalues == dnorm(1)), labels=c(-1, 1)) axis(1, at=which(dvalues == dnorm(2)), labels=c(-2, 2)) ``` As you can see, `dnorm` will give us the "height" of the pdf of the normal distribution at whatever Z-score we provide as an argument to `dnorm`. --- ## `pnorm` The function `pnorm` returns the integral from $-\infty$ to $q$ of the pdf of the normal distribution where $q$ is a Z-score. Try to guess the value of `pnorm(0)`. (`pnorm` has the same default `mean` and `sd` arguments as `dnorm`). ```{r} # To be clear about the arguments in this example: # q = 0, mean = 0, sd = 1 pnorm(0) ``` The `pnorm` function also takes the argument `lower.tail`. If `lower.tail` is set equal to `FALSE` then `pnorm` returns the integral from $q$ to $\infty$ of the pdf of the normal distribution. Note that `pnorm(q)` is the same as `1-pnorm(q, lower.tail = FALSE)` ```{r} pnorm(2) pnorm(2, mean = 5, sd = 3) pnorm(2, mean = 5, sd = 3, lower.tail = FALSE) 1 - pnorm(2, mean = 5, sd = 3, lower.tail = FALSE) ``` `pnorm` is the function that replaces the table of probabilites and Z-scores at the back of the statistics textbook. Let's take our vector of Z-scores from before (`z_scores`) and compute a new vector of "probability masses" using `pnorm`. Any guesses about what this plot will look like? ```{r} pvalues <- pnorm(z_scores) # Now we'll plot these values plot(pvalues, # Plot where y = values and x = index of the value in the vector xaxt = "n", # Don't label the x-axis type = "l", # Make it a line plot main = "cdf of the Standard Normal", xlab= "Quantiles", ylab="Probability Density") # These commands label the x-axis axis(1, at=which(pvalues == pnorm(-2)), labels=round(pnorm(-2), 2)) axis(1, at=which(pvalues == pnorm(-1)), labels=round(pnorm(-1), 2)) axis(1, at=which(pvalues == pnorm(0)), labels=c(.5)) axis(1, at=which(pvalues == pnorm(1)), labels=round(pnorm(1), 2)) axis(1, at=which(pvalues == pnorm(2)), labels=round(pnorm(2), 2)) ``` It's the plot of the cumulative distribution function of the normal distribution! Isn't that neat? --- ## `qnorm` The `qnorm` function is simply the inverse of the cdf, which you can also think of as the inverse of `pnorm`! You can use `qnorm` to determine the answer to the question: What is the Z-score of the $pth$ quantile of the normal distribution? ```{r} # What is the Z-score of the 50th quantile of the normal distribution? qnorm(.5) # What is the Z-score of the 96th quantile of the normal distribution? qnorm(.96) # What is the Z-score of the 99th quantile of the normal distribution? qnorm(.99) # They're truly inverses! pnorm(qnorm(0)) qnorm(pnorm(0)) ``` Let's plot `qnorm` and `pnorm` next to each other to further illustrate the fact they they are inverses. ```{r, warning=FALSE} # This is for getting two graphs next to each other oldpar <- par() par(mfrow=c(1,2)) # Let's make a vector of quantiles: from 0 to 1 by increments of .05 quantiles <- seq(0, 1, by = .05) quantiles # Now we'll find the Z-score at each quantile qvalues <- qnorm(quantiles) qvalues # Plot the z_scores plot(qvalues, type = "l", # We want a line graph xaxt = "n", # No x-axis xlab="Probability Density", ylab="Z-scores") # Same pnorm plot from before plot(pvalues, # Plot where y = values and x = index of the value in the vector xaxt = "n", # Don't label the x-axis type = "l", # Make it a line plot main = "cdf of the Standard Normal", xlab= "Quantiles", ylab="Probability Density") # These commands label the x-axis axis(1, at=which(pvalues == pnorm(-2)), labels=round(pnorm(-2), 2)) axis(1, at=which(pvalues == pnorm(-1)), labels=round(pnorm(-1), 2)) axis(1, at=which(pvalues == pnorm(0)), labels=c(.5)) axis(1, at=which(pvalues == pnorm(1)), labels=round(pnorm(1), 2)) axis(1, at=which(pvalues == pnorm(2)), labels=round(pnorm(2), 2)) # Restore old plotting settings par(oldpar) ``` --- ## `rnorm` If you want to generate a vector of normally distributed random numbers, `rnorm` is the function you should use. The first argument `n` is the number of numbers you want to generate, followed by the standard `mean` and `sd` arguments. Let's illustrate the weak law of large numbers using `rnorm`. ```{r, warning=FALSE} # set.seed is a function that takes a number as an argument and sets a seed from # which random numbers are generated. It's important to set a seed so that your # code is reproduceable. If you wanted to you could always set your seed to the # same number. I like to set seeds to the "date" which is really just # the arithmetic equation "month minus day minus year". So today's seed # is -2006. set.seed(10-1-2015) rnorm(5) # If I set the seed to the same seed again, I'll generate the same vector of # numbers. set.seed(10-1-2015) rnorm(5) # Now onto using rnorm # Let's generate three different vectors of random numbers from a normal # distribution n10 <- rnorm(10, mean = 70, sd = 5) n100 <- rnorm(100, mean = 70, sd = 5) n10000 <- rnorm(10000, mean = 70, sd = 5) # Let's just look at one of the vectors n10 ``` Which historgram do you think will be most centered around the true mean of 70? ```{r, warning=FALSE} # This is for getting two graphs next to each other oldpar <- par() par(mfrow=c(1,3)) # The breaks argument specifies how many bars are in the histogram hist(n10, breaks = 5) hist(n100, breaks = 20) hist(n10000, breaks = 100) # Restore old plotting settings par(oldpar) ``` --- ## Closing thoughts These concepts generally hold true for all the distribution functions built into R. You can learn more about all of the distribution functions by typing `help(Distributions)` into the R console. If you have any questions about this demonstration or about R programming please send me an [email](http://seankross.com/about/). If you'd like to change or contribute to this document I welcome pull requests on [GitHub](https://github.com/seankross/seankross.github.io/tree/master/notes). This document and all code contained within is licensed [CC0](https://creativecommons.org/publicdomain/zero/1.0/).