---
title: "Math 189/289 Lab 4"
author: "Jiaqi Guo"
output:
  html_notebook: default
  pdf_document: default
---

# Create Data
Let's say today we are going to work with some simulated data. We generate $N$ numbers from bernoulli distribution with success rate 0.3, and take them as our population.
```{r}
N <- 1000
data.population <- rbinom(n=N, size=1, prob=0.3)
```

One can give a story line to this data. If put in the same context as the video games dataset that we have for homework, for example, this could be the response of the students in the whole school whether they played games in the past week or not. Thus, the school has 1000 students in total. 1 indicates the student played, and 0 indicates did not play.

Then in that sense, we might not have observed every one of these responses as in for the survey data in the video games dataset. Thus, we sample, say $n$, observations from the population, and take them as our sample units.
```{r}
n <- 300
ind.sample <- sample.int(n=N, size=n)
data.sample <- data.population[ind.sample]
```

# Sample Statistics
Now we will stick with the story line. Once we have the sample data, we can compute a point estimate, as well as an interval estimate for the fraction of the students who played games in the past week or not.

A point estimate in this case is just the mean of the sample,
```{r}
mean.sample <- mean(data.sample)
mean.sample
```

To get an interval estimate for the fraction, we follow the derivation in the lecture slides, the interval estimate is then given by,
$$ \left(\bar x - 1.96 \sqrt{\frac{\bar x (1 - \bar x)}{n-1} \frac{N - n}{N}}, \bar x + 1.96 \sqrt{\frac{\bar x (1 - \bar x)}{n-1} \frac{N - n}{N}} \right), $$
where $\bar x$ indicates the sample mean, $N$ indicates the population size, and $n$ indicates the sample size. Thus, an interval estimate is then given as the following.
```{r}
width <- 1.96 * sqrt(mean.sample*(1-mean.sample)*(N-n)/((n-1)*N))
int.sample <- c(mean.sample - width, mean.sample + width)
int.sample
```

# Bootstrap Estimate
Another popular method is using the idea of bootstrap. We need a bootstrap population to start with. With the population size of $N=1000$, and sample size of $n=300$, approximately each sample occurs about 3 times in the bootstrap population. One can simply duplicate each observation in the sample 3 times, and treat the resulting sample as the bootstrap population.

One can also sample from the sample with replacement.
```{r}
ind.boot <- sample.int(n, size=N, replace=TRUE)
data.boot <- data.sample[ind.boot]
```

With the bootstrap population, we are now ready to generate bootstrap sample means. We will take, say 2000, random samples of size $n$.
```{r}
B <- 2000
boot.sample.mean <- rep(NA, B)
for(i in 1:B){
  ind <- sample.int(N, size=n, replace=FALSE)
  boot <- data.boot[ind]
  boot.sample.mean[i] <- mean(boot)
}
```

Let's take a look at the distribution of bootstrap sample means. 
```{r}
hist(boot.sample.mean)
```

The point estimate is the mean of bootstrap sample means.
```{r}
mean.boot <- mean(boot.sample.mean)
mean.boot
```

There are two approaches one can take from here. We can derive a confidence interval using the bootstrap sample means, as given by
$$ \left( \bar x - 1.96 s, \bar x + 1.96s \right) $$
```{r}
s <- sd(boot.sample.mean)
int.boot <- c(mean.boot - 1.96*s, mean.boot + 1.96*s)
int.boot
```

Another approach to derive an interval estimate using the bootstrap sample means, one can simply extract the 0.025-quantile and 0.975 quantile of the bootstrap sample means and arrive at an interval estimate.
```{r}
int.boot <- c(quantile(boot.sample.mean, 0.025), quantile(boot.sample.mean, 0.975))
int.boot
```

# Comparing Two Distributions
To compare two distributions, the Kolmogorov-Smirnov (KS) statistics is a helpful measure. For example, say if we take another sample from the population data earlier, and compare it with the earlier sample we took.
```{r}
ind.sample2 <- sample.int(n=N, size=n)
data.sample2 <- data.population[ind.sample2]
ks.test(data.sample, data.sample2)
```

Now let's take another sample that is completely different.
```{r}
data.sample2 <- rbinom(n=n, size=1, prob=0.7)
ks.test(data.sample, data.sample2)
```

# Important Features
Often we woule like to pick out important reasons resulting an outcome. For example, we would like to figure out what are the import reasons why students like/dislike video games.

We first load the data.
```{r}
data <- read.table("videodata.txt", header=TRUE)
```

The package we are going to introduce today is named 'tree'.
```{r}
install.packages("tree")
library(tree)
```

Since our task only ask for like or dislike, we need to congregate the column 'like' into 2 categories from 5, as a new column named 'dis_like'. We will denote 1=never played, 4=not really, and 5=not at all of the column 'like' as 0=dislike in the new column 'dis/like', and the rest as 1=like.
```{r}
data['dis_like'] <- rep(NA, dim(data)[1])
for(i in 1:dim(data)[1]){
  like <- data[i, 'like']
  if(like==0 || like==4 || like==5){
    data[i, 'dis_like'] = 0
  }else{
    data[i, 'dis_like'] = 1
  }
}
```

We then fit a tree on the data frame with appropriate columns.
```{r}
data.tree <- tree(dis_like~educ+sex+age+home+math+work+own+cdrom+grade, data=data)
plot(data.tree, type="uniform")
text(data.tree)
```

