--- title: "lab5_start_code" author: "Kai Ping (Brian) Leung" date: "10/25/2019" output: pdf_document editor_options: chunk_output_type: console --- {r setup, include=FALSE} knitr::opts_chunk\$set(echo = TRUE)  # Prerequisite {r} library(tidyverse) theme_set(theme_classic())  # Homework 2 Recap $$\begin{split} \Pr (\textbf{y} | \pi) & = \prod_{i = 1}^{n} \pi^{y_i} (1-\pi)^{1-y_i} \\ L(\pi | \textbf{y}) & \propto \prod_{i = 1}^{n} \pi^{y_i} (1-\pi)^{1-y_i} \\ & \propto \pi^{y_1} \times \pi^{y_2} \times \dots \times \pi^{y_n} \times (1-\pi)^{1- y_1}\times (1-\pi)^{1-y_2}\times \dots \times (1-\pi)^{1-y_{n}} \\ & \propto \pi^{y_1 + y2 + \dots + y_n} \times (1-\pi)^{n - (y_1 + y2 + \dots + y_n)} \\ & \propto \pi^{\mathbf{y}} \times (1-\pi)^{n - \mathbf{y}} \\ \\ \text{Given the actual data are } &(1, 0, 0, 1, 0, 0, 0, 0), \text{ i.e.,} \\ \mathbf{y} &= 2 \hskip{0.5in} n - \mathbf{y} = 8-2 = 6 \\ \\ \therefore L(\pi | \textbf{y}) & \propto \pi^{2} \times (1-\pi)^{6} \end{split}$$ # Implement the likelihood estimator in R {r} bern_table <- tibble(pi = seq(from = 0, to = 1, by = 0.01), llk = pi ^ 2 * (1-pi) ^ 6) ggplot(bern_table, aes(x = pi, y = llk))+ geom_line()+ labs(title = "Likelihood function of pi") logit.llk <- function(beta, x , y) { x <- as.matrix(x) xb <- x %*% beta pi <- 1 / (1 + exp(-xb)) prod(pi ^ y * (1 - pi) ^ (1 - y)) } y <- c(1, 0,0,0,1,0,0,0) x <- c(5, -2, -2,-2,5,-2,-2,-2) betas <- -10:10 bern.table <- tibble( beta = betas, llk = map_dbl(betas, logit.llk, y = y, x = x) ) ggplot(bern.table, aes(x = beta, y = llk)) + geom_line()  # What if 5 success out of 8? {r} bern_table2 <- tibble(pi = seq(from = 0, to = 1, by = 0.01), llk = pi ^ 5 * (1-pi) ^ 3) ggplot(bern_table2, aes(x = pi, y = llk))+ geom_line()+ labs(title = "Likelihood function of pi")  # What if 7 success out of 8? {r} bern_table3 <- tibble(pi = seq(from = 0, to = 1, by = 0.01), llk = pi ^ 7 * (1-pi) ^ 1) ggplot(bern_table3, aes(x = pi, y = llk))+ geom_line()+ labs(title = "Likelihood function of pi")  # Combine {r} bern_all <- bind_rows( bern_table %>% mutate(success = "2 success"), bern_table2 %>% mutate(success = "5 success"), bern_table3 %>% mutate(success = "7 success") ) ggplot(bern_all, aes(x = pi, y = llk))+ geom_line() + labs(title = "Likelihood function of pi") + facet_wrap(~ success, scale = "free_y")