---
title: 'CSSS/POLS 510 MLE'
subtitle: "Lab 3. Heteroskedastic Normal"
author: Minji Jeong
date: \today
classoption: xcolor=dvipsnames
output:
  beamer_presentation:
    colortheme: seahorse
    fonttheme: structurebold
    theme: CambridgeUS
    highlight: default
linkcolor: blue
editor_options: 
  chunk_output_type: console
header-includes:
 \usepackage{float}
 \usepackage{xcolor}
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,message = FALSE,warning = FALSE,error = FALSE)
library(tidyverse);library(MASS);library(simcf)
```


# Housekeeping

- All Problem Set answers should be submitted on Canvas. 

    - Submit your file in pdf form and make sure you present the R code you wrote (echo=TRUE). 
    - Project proposal, final paper, and poster should be sent to Chris's email.

- PS 1 grade will be released next week. 

- For office hours, please contact me 24 hours before the time you want to have a zoom meeting. 

- Make sure you install Chris's [simcf](https://faculty.washington.edu/cadolph/index.php?page=60) package. 


# Least Squares

::: columns
:::: column
\scriptsize
- Linear homoskedastic:
\begin{align*}
Y_i = x_i \mathcal{\beta}  + \varepsilon, \quad  \varepsilon \sim N(0, \sigma^2)
\end{align*}
- Estimating the slope:
\begin{align*}
\quad \hat{\beta}_j &= \frac{\text{Cov}(X_j, Y)}{\text{Var}(X_j)}
\end{align*}
- Matrix Algebra:
\begin{align*}
\quad \hat{\beta} &= (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}
\end{align*}
::::
:::: column
```{r, echo=F, fig.width=6 , fig.height=6/1.618}
swiss %>%
  ggplot(aes(x = Education, y = Fertility)) +
  geom_point(size=2) +
  geom_smooth(method = "lm", 
              formula = "y~x",
              se = FALSE, 
              color = "blue") +
  labs(y = "Fertility", x = "% Education Beyond Primary School") +
  theme_minimal() +
  theme(axis.text = element_text(size = 15, color = "black"))
```
::::
:::


# Least Squares

::: columns
:::: column

\scriptsize

\begin{align}
\hat{y_i} & = \hat{\beta_0} + \hat{\beta_1} x_i \\
\hat{e_i} & = y_i - \hat{y}_i \\
\sum \hat{e_i}\textcolor{red}{^2} & = \sum (y_i - \hat{y}_i)\textcolor{red}{^2}\\
\sum \hat{e_i}\textcolor{red}{^2} & = \sum (y_i - \hat{\beta_0} - \hat{\beta_1} x_i)\textcolor{red}{^2} 
\end{align}

\tiny

Note: (3) and (4) are equivalent, and provide the *Residual Sum of Squares* (RSS) or *Sum of Squared Residuals* (SSR).

\vspace{5mm}

Choosing the best combination of $\hat{\beta_j}$ that minimizes the SSR provide the best fit for the model.

::::
:::: column
```{r, echo=F, fig.width=6 , fig.height=6/1.618}
swiss %>%
  lm(Fertility ~ Education, data = .) %>%
  broom::augment() %>%
  ggplot(aes(x = Education, y = Fertility)) +
  geom_point(size=2) +
  geom_smooth(method = "lm", 
              formula = "y~x",
              se = FALSE, 
              color = "blue") +
  labs(y = "Fertility", x = "% Education Beyond Primary School") +
  geom_segment(aes(yend = .fitted, xend = Education), linetype = 3,size=1,color="red") +
  theme_minimal() +
  theme(axis.text = element_text(size = 15, color = "black"))
```
::::
:::


# Maximum Likelihood Estimation

- How do we estimate the MLE?
  1. Define a probability model (PDF): $Y_i \sim N(\mu_i, \sigma^2)$.
  2. Derive the log-likelihood function.
  3. Reduce to sufficient statistics and substitute systematic component.
  4. Use `optim()` or any other function to find the maxima.



# Normal homoskedastic

Two **different** notations for the **same** model.

\vspace{5mm}

::: columns
:::: column
LS notation:
\scriptsize
$$
\begin{aligned}
\varepsilon_i & \sim N(0, \sigma^2) & (stochastic) \\
Y_i & = x_i \mathcal{\beta} & (systematic) \\
Y_i & = x_i \mathcal{\beta}  + \varepsilon  & (stochastic + systematic) \\
\end{aligned}
$$


::::
:::: column
MLE notation:
\scriptsize
$$
\begin{aligned}
Y_i & \sim N(\mu_i, \sigma^2) & (stochastic) \\
\mu_i & = x_i \beta & (systematic) \\
Y_i & \sim N( x_i \beta, \sigma^2)  & (stochastic + systematic) \\
\end{aligned}
$$
::::
:::





# MLE general notation

Systematic–Stochastic Decomposition:



$$
\begin{aligned}
Y_i & \sim f(\theta_{i}, \alpha)  & \text{(stochastic)}\\
\theta_{i} & = g(\boldsymbol{x}_{i} \boldsymbol{\beta}) & \text{(systematic})
\end{aligned}
$$
where

  - $Y_i$ is a outcome random variable.
  - $f()$ is a probability density function.
  - $\theta_i$ is a systematic feature of the PDF that varies over i.
  - $\alpha$ is an ancillary parameter (feature of f that we treat as constant).
  - $g()$ functional form for reparametrization of the data model.
  - $x_i$ explanatory variables vector.
  - $\beta$ vector of effect parameters.




# MLE - Homoskedastic normal


```{r,echo=FALSE, fig.width=4.5 , fig.height=4.5/1.618, fig.align='center'}

## set simulation parameters
mu = 67
sigma2 = 10
n = 100

## generate random normally distributed samples
y = rnorm(n, mu, sd = sqrt(sigma2))

## Hypothetical values of mu
mu_range <- seq(from=55, 
                to=75,
                length.out= 1000)

# Let's visualize the profile likelihoods and the SSR


## function that will iterate SSR for each value of mu

SSR_profile <- function(mu, y) {
  # create empty vector
  SSR.vector <- vector(mode = "numeric", length = length(mu))
  
  for (i in 1:length(SSR.vector)){
    # estimate SSR
    SSR.vector[i] <- sum((y - mu[i]) ^ 2)
  }
  
  return(SSR.vector)
}

## Use optim to find the minima
optim_mu <- optim(par=10,
                  fn=SSR_profile,
                  method = "BFGS",
                  y=y)


## function for log likelihood

loglike_norm <- function(mu, sigma2, y) {
  # create empty vector
  likelihood.vector <- vector(mode = "numeric", length = length(mu))
  
  # for each mu, calculate the likelihood value
  for (i in 1:length(likelihood.vector)) {
    likelihood.vector[i] <- sum(dnorm(x = y, 
                                      mean = mu[i], 
                                      sd = sqrt(sigma2), 
                                      log = T))
  }
  
  # return vector of likelihood values spanning range of mu
  return(likelihood.vector)
}

df_vis <-
  tibble(SSR = SSR_profile(mu=mu_range,y=y),
         LogLikelihood_Normal =loglike_norm(mu=mu_range, sigma2 = sigma2, y = y),
         mu = mu_range)

df_vis %>% 
  rename("SSR (LS)" = SSR,
         "Log-Likelihood Normal (MLE)"=LogLikelihood_Normal) %>% 
  pivot_longer(cols = c("SSR (LS)","Log-Likelihood Normal (MLE)"),
               names_to = "profile",
               values_to = "value") %>% 
  ggplot(aes(y=value,
             x=mu)) +
  geom_line() +
  geom_vline(xintercept = optim_mu$par,
             linetype = "dashed",
             color="red") +
  facet_wrap(~profile, scales = "free") +
  theme_minimal()+
  theme(
    panel.background = element_rect(fill = NA),
    axis.line.x.bottom = element_line(size = 0.5),
    axis.ticks.length = unit(0.5, "char"),
    legend.position = "none",
    axis.text = element_text(size = 7, color = "black")
  )

```



# 2. Heteroskedastic normal

+ Steps
  + The full `R` code can be found [\underline{here from Chris' website}](https://faculty.washington.edu/cadolph/mle/HeteroskedMLE.r)

    - Generate Data

    - Fit OLS  - `lm()`

    - Fit MLE  - `optim()`






# 

\Huge

\centerline{F I N}

