Goal: After completing this lab, you should be able to…

• Fit any well-known distribution (that is implemented in R) using either maximum likelihood or method of moments estimators.
• Plot a fitted distribution by adding it to a histogram of data that was modeled.
• Check how well a distribution models some data using a QQ-Plot.

In this lab we will use, but not focus on…

• R Markdown. This document will serve as a template. It is pre-formatted and already contains chunks that you need to complete.
• Importing data. Code to import the necessary data is provided.

Some additional notes:

• Please see Carmen for information about submission, and grading.
• You may use this document as a template. You do not need to remove directions. Chunks that require your input have a comment indicating to do so.
• The reading posted with the lab and the associated R Markdown file contain details on performing the required tasks through a number of examples.
• Recall the tutorial information posted in Lab 01. It may be extremely useful here.

# Exercise 1 - Sleep Data, Beta Distribution

For this exercise, we will return to the sleep data from class. Please refer to the notes for a complete description of this dataset.

library(readr)
fitbit_sleep = read_csv("https://daviddalpiaz.github.io/stat3202-au18/data/fitbit-sleep.csv")

The above chunk will read the data as a data frame named fitbit_sleep. (Actually, a tibble, but that is unimportant at the moment.) In general, you should avoid using an absolute reference to data on your own machine. Instead, you should use a relative reference, and include the data with your R Markdown document. Here we avoid the issue altogether and simply read the data directly from the web. This way, the above code will work on anyone’s machine.

percent_REM = fitbit_sleep$min_REM / fitbit_sleep$min_sleep

For this exercise, we will investigate a new variable called percent_REM that is created above. This variable measures the proportion of time asleep spent in the REM sleep stage.

We will consider two probability models for this variable, a beta model and a normal model.

range(percent_REM)
## [1] 0.07474747 0.28508772

Notice that the range of the observed data that is calculated above. Furthermore, since this variable is a proportion, it could take values between 0 and 1. (Although, it is rather unlikely that you would spend all night in REM sleep.) If we think of the sample space of the continuous distributions that we know, we should remember that the beta distribution takes as input values between 0 and 1. This makes beta a good candidate.

Recall, that a random variable $$X$$ that follows a beta distribution has probability density function

$f(x \mid \alpha, \beta) = \left[\frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)}\right] x^{\alpha - 1}(1 - y)^{\beta - 1}, 0 \leq x \leq 1, \alpha > 0, \beta > 0$

Also, the beta distribution has

$\text{E}[X] = \frac{\alpha}{\alpha + \beta}$

and

$\text{Var}[X] = \frac{\alpha\beta}{(\alpha + \beta)^2(\alpha + \beta + 1)}$

While it is not possible to find an analytic solution for the MLE, given $$X_1, X_2, \ldots X_n$$ assumed to be a random sample from a beta distribution, we can use method of moments to obtain estimators:

$\tilde{\alpha} = \bar{x}\left(\frac{\bar{x}(1 - \bar{x})}{s^2} - 1\right)$

$\tilde{\beta} = (1 - \bar{x})\left(\frac{\bar{x}(1 - \bar{x})}{s^2} - 1\right)$

where

$\bar{x} = \frac{1}{n}\sum_{i = 1}^{n}x_i$

$s^2 = \frac{1}{n}\sum_{i = 1}^{n}(x_i - \bar{x})^2$

Do the following:

• Calculate estimates for the parameters assuming a beta distribution using this data. (Use the provided method of moments estimators above.)
• Add the fitted model to the provided histogram. (That is, add the density for the specific beta distribution you found.)
• Create a QQ-plot to better determine how well this distribution fits the data.
# perform estimate calculations here
# remove these comments and insert your code in this chunk
hist(percent_REM, xlim = c(0, 1), ylim = c(0, 9), probability = TRUE,
main = "Histogram of Percent REM Sleep",
xlab = "Percent REM Sleep")
box()
grid()