---
title: "Math 189/289 Lab 3"
author: "Alex Hanbo Li"
output:
  html_notebook: default
  pdf_document: default
---

# Dataframe
R has a working directory. This should be the main location to store external data files, as well as other code files.
```{r}
getwd()
```
One can also set the working directory to a desired location.
```{r}
setwd() # desired directory as argument
getwd()
```
If you are not sure where the source file is located, you can use "Session --> Set Working Directory --> To Source File Location".

First, let's load in the video game data for homework 2.
```{r}
setwd("~/Dropbox/Math189Winter2017/lab3") # change to your own working directory
data <- read.table("videodata.txt", header=TRUE)
head(data)
```

If a question was not answered or improperly answered, it was decoded as 99. There are total number of 65 such fields.
```{r}
data[data == 99] <- NA
sum(is.na(data))
```

## `attach()` and `detach()` functions
The **$** notation for list components is not always very convenient.
Sometimes it is more convenient to temporarily make the components of data frame
visible as variables under their component name. We can use `attach()` function in R to do this.

For example, if we want to access the **age** variable, before attaching the data frame `data`, we can only access by:
```{r}
data$age # returns an array
data['age'] # return a dataframe
```

If we directly type in `age`, an error will be returned.
```{r}
age
```

After attaching the dataframe, we have an individual variable called `age`.
```{r}
attach(data)
age
```

However, now any change on `age` will only affect the new variable `age` in the working directory, not the `age` attribute of the data frame `data`. 
```{r}
age[1] <- 1000
age
data$age
```

To detach a data frame, use function `detach()`.
```{r}
search()
ls(2)
detach(data)
search()
```

## More on estimating the distribution of data
### Histogram
```{r}
summary(data$age)
hist(data$age)
hist(data$age, breaks = 25, probability = TRUE, density = 20, col = 3, border = 3)
lines(density(data$age, bw = 0.25), col = 1)
lines(density(data$age, bw = 0.5), col = 2)
rug(data$age)
```

### Empirical cumulative distribution function (`ecdf()`)
```{r}
plot(ecdf(data$age), do.points = F, verticals = T)
x <- seq(min(data$age)-4, max(data$age), length.out = 1000)
lines(x, pnorm(x, mean=mean(data$age), sd=sd(data$age)), lty=3, col = 2)
```
Let's compare it with a Normal distribution with the same sample mean and variance (red dotted line above).

Of course, one would check the Normality more carefully using Q-Q plot. We see a very long right tail.
```{r}
par(pty = 's') # plot in a square
qqnorm(data$age)
qqline(data$age)
```

More formally, we can do a Shapiro-Wilk test. In this case, the p-value is significantly small and we can reject the Normality hypothesis $H_0$.
```{r}
shapiro.test(data$age)
```

# Sampling and Bootstrap
### Sampling
Sampling (with or without replacement) can be done in R using the function `sample()`.
```{r}
sample(1:10, size = 10, replace = F)
sample(1:10, size = 10, replace = T)
```

### Bootstrap
Suppose we want to find the point estimate of the proportion of male in the population and its confidence interval. Note that the *sex* variable is 1 if the student is male and 0 if female. Then the point estimation of the proportion of male is:
```{r}
male.percentage <- mean(data$sex)
male.percentage
```

Now we also want to have a confidence interval of this estimator. However, clearly the distribution of *sex* variable is not Normal, it is a Bernoulli random variable. We know our data were drawn from a population with size $N = 314$. Hence, we first create a bootstrap population of this size by repeating every sample for $\frac{314}{91} = 3.45$ times. Here, we'll just specify the parameter *length.out* to be 314.
```{r}
boot.population <- rep(data$sex, length.out = 314)
length(boot.population)
```

Then we will choose $n = 91$ samples from the Bootstrap population and call this a Bootstrap sample.
```{r}
sample1 <- sample(boot.population, size = 91, replace = FALSE)
```

Continue this procedure until we have 400 Bootstrap samples.
```{r}
set.seed(189289)
B = 400 # the number of bootstrap samples we want
boot.sample <- array(dim = c(B, 91))
for (i in 1:B) {
  boot.sample[i, ] <- sample(boot.population, size = 91, replace = FALSE)
}
```

Then we can calculate the sample mean for each Bootstrap sample (i.e. each row of the Bootstrap sample matrix).
```{r}
boot.mean <- apply(X = boot.sample, MARGIN = 1, FUN = mean)
head(boot.mean)
```

Let's see the histogram of these Bootstrap sample means.
```{r}
hist(boot.mean, breaks = 20, probability = TRUE, density = 20, col = 3, border = 3)
lines(density(boot.mean, adjust = 2), col = 2)
```

Check Normality by Q-Q plot and Shapiro-Wilk test.
```{r}
par(pty = 's')
qqnorm(boot.mean)
qqline(boot.mean)

shapiro.test(boot.mean)
```

So we can accept that the sample mean follows a Normal distribution. Then we can construct $95\%$ confidence intervals.
```{r}
boot.sd <- sd(boot.mean)
male.percentage + c(-1, 1)*1.96*boot.sd
```
