---
title: "Lab 2: Practice Code"
subtitle: "Probability Distributions, Statistical Inference, and Ordinay Least Squares"
author: "Minji Jeong"
date: \today
output: pdf_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Prerequisite
```{r, warning=F, message=FALSE}
# Good practice to remove all objects from the workspace
rm(list = ls())

# Use library() for packages you need, or source() for other R files.
library(tidyverse)

# Setting the seed ensures that we get the same random draw over and over again.
set.seed(2025)

rnorm(5) # Check
```



# Probability Distributions

## 0. Calculate the following operations by hand (... meaning by R)
a) 
$$ \sum\limits_{i=1}^{5} i = $$
b) 
$$ \prod\limits_{i=1}^{5} i = $$
c) 
$$ 5! \times 10^{3!} \times e^{4} = $$

```{r}

# a) 
sum(1:5)

# b)

prod(1:5)
#Check 
1*2*3*4*5

# c)

#factorial(5)
#5*4*3*2*1

factorial(5) * 10^factorial(3) * exp(4)


```



## 1. Build a Bernoulli distribution using the `sample()` function, where the probability of *success* is 0.7. Run `?sample` if you are unsure how the function works.

```{r}
# Create an imaginary person to flip the coin once for you 
sample(x = c(0, 1),
       size = 1,
       prob = c(0.3, 0.7)
       )


```

## 2. How do you know if it is working properly? Conduct simulation with `for loop` to check if the assigned probabilities are matched with the empirics 
```{r}
# Specify the number of simulations
sims <- 10000

# Specify the probability
ProbSuccess <- 0.7

# Create an empty vector as "container"
BernResult <- vector(mode = "numeric", 
                     length = sims)

# For loop
for (i in 1:sims) {
  BernResult[i] <- sample(c(0, 1), 
                           size = 1, 
                           prob = c(1 - ProbSuccess,
                                    ProbSuccess)
                           )
}

mean(BernResult)

```

```{r, results=FALSE}
## Faster way w/o loop...but we want to build the intuition
 sample(c(0, 1),
        size = sims,
        replace = TRUE,
        prob = c(0.3, 0.7)
        )
## Or use rbinom (rmb: bernoulli is just a special case of binomial when N = 1)
 rbinom(sims,
        size = 1,
        prob = 0.7)

```

## 3. Plot the above Bernoulli distribution
```{r, warning=F, message=F, fig.width=5, fig.height=5/1.6, fig.align='center'}
# Base graphics
hist(BernResult)

# ggplot2
BernDF <- tibble(Outcome = BernResult) 


ggplot(data=BernDF, 
       aes(x=Outcome)) +
  geom_histogram() + 
  scale_x_continuous(breaks = c(0, 1), expand = c(1, 0)) +
  labs(y = "Frequency", 
       x = "x",
       subtitle = "Bernoulli Distribution: Frequency")+
  theme_bw()


# Adding: aes(y = stat(count / sum(count))) in geom_hist ; What does it do?
ggplot(data=BernDF, 
       aes(x=Outcome)) +
  geom_histogram(aes(y = after_stat(count/sum(count)))) +   
  scale_x_continuous(breaks = c(0, 1), expand = c(1, 0)) +
  labs(y = "Prob", x = "x",
       subtitle = "Bernoulli Distribution: Probability")+
  theme_bw()
 

# y = stat() means...
BernResult |> 
  as_tibble() |> 
  count(value) |> 
  mutate(TotalNum = sum(n),
         Prob = n/TotalNum)


```

## 4. Based on the above, generate a binomial distribution, with number of trials equal to 10, without using rbinom()
```{r}
# Create an imaginary person to flip the coin ten times for you
# Let's test it outside of the loop: 
sample(0:1,
       size = 10,
       replace = TRUE,
       prob = c(1- ProbSuccess,
                ProbSuccess)
       )

# Create number of simulations and an empty vector as container 
BinoResult <- vector(mode = "numeric", length = sims)

for (i in 1:sims) {
  
  # Create an imaginary person to flip the coin ten times for you
  flips <- sample(c(0, 1),
                  size = 10,
                  replace = TRUE, 
                  prob = c(1- ProbSuccess, 
                           ProbSuccess)
                  )
  
  # Sum up the number of "success" for that person
  count <- sum(flips)
  
  # Store it into the container; repeat 10,000
  BinoResult[i] <- count
  
}

# Faster way w/o loop and sample()...but we want to build the intuition
# rbinom(n = sims,
#        size = 10,
#        prob = ProbSuccess)
```

## 5. Plot the above binomial distribution
```{r}
# Base graphics
hist(BinoResult)

# ggplot2
BinoResult |> 
  as_tibble() |>  # Convert to data frame
  ggplot(aes(x = value)) +
  geom_histogram(aes(y = stat(count/sum(count)))) +   
  scale_x_continuous(breaks = 1:10) +
  labs(y = "Prob", x = "x",
       subtitle = "Binomial Distribution")+
  theme_bw()

```


```{r, include=F}
 rbinom(n = sims,
        size = 10,
        prob = ProbSuccess) |> 
  as_tibble() |>  # Convert to data frame
  ggplot(aes(x = value)) +
  geom_histogram(aes(y = stat(count/sum(count)))) + 
  scale_x_continuous(breaks = 1:10) +
  labs(y = "Probability", x = "x",
       subtitle = "Binomial Distribution")+
  theme_bw()


```

## 6. Explore the rbinom, dbinom, pbinom functions. What do they do? Answer the following questions:
  a) The probability of a coin landing on head is 0.7. If you were to flip the coin 10 times, what is the probability of getting exactly 7 heads?
  b) What is the probability of getting 7 heads or less?
  c) How do you know (b) is true?
```{r}
# a) Pr(exactly 7 heads) -> PDF
dbinom(x = 7, size = 10, prob = 0.7)
  
# b) Pr(7 heads or less) -> CDF
pbinom(q = 7, size = 10, prob = 0.7)

pbinom(q = 7, size = 10, prob = 0.7, lower.tail = FALSE)

1 - pbinom(q = 7, size = 10, prob = 0.7)


# c) Double check
dbinom(x = 0:7, size = 10, prob = 0.7) |> 
  sum()

```



-----------


## 4. Based on the above, generate a binomial distribution, with number of trials equal to 10, without using rbinom()
```{r}
# Imagine you flip the coin ten times
# Let's test it outside of the loop: 

sample(x = c(0,1),
       size = 10,
       replace = T,
       prob = c(1-ProbSuccess,ProbSuccess))
```

```{r}
# Create number of simulations and an empty vector as container 
set.seed(12345)

sims <- 9000

BinoResult <- vector(mode = "numeric", length = sims)

for (i in 1:sims) {
  
  # Flip the coin ten times
  flips <- sample(x = c(0,1),
                  size = 10,
                  replace = T,
                  prob = c(1-ProbSuccess,ProbSuccess))

  # Store the sum of the number of "success" into the container; repeat 9000
  BinoResult[i] <- sum(flips)
}


head(BinoResult)
 # What are we doing with the following code?
BinoResult7 <- sum(BinoResult == 7)
BinoResult7/sims
```


```{r}
## EXTRA

# Create number of simulations  
sims <- 9000

# Get 9000 binomial draws (each draw = #successes in 30 trials)
BinoResult <- rbinom(n = sims, size = 10, prob = ProbSuccess)

BinoResult7 <- sum(BinoResult == 7)
BinoResult7/sims
```


## 5. Plot the above binomial distribution
```{r}

hist(BinoResult)

```

## 6. Explore the rbinom, dbinom, pbinom functions. What do they do? Answer the following questions:
  a) The probability of a coin landing on head is 0.7. If you were to flip the coin 10 times, what is the probability of getting exactly 7 heads?
  b) What is the probability of getting 7 heads or less?
  c) How do you know (b) is true?
```{r}
# a) Pr(exactly 7 heads) -> PdF
dbinom(x = 7, size = 10, prob = 0.7)
  
# b) Pr(7 heads or less) -> CDF
pbinom(q = 7, size = 10, prob = 0.7)

1-pbinom(q = 7, size = 10, prob = 0.7)
pbinom(q = 7, size = 10, prob = 0.7, lower.tail = F)

# c) Double check
dbinom(x = 0:7, size = 10, prob = 0.7) |> 
  sum()
```












