% \documentclass[serif]{beamer} % Serif for Computer Modern math font.
\documentclass[serif, handout]{beamer} % Handout mode to ignore pause statements
\hypersetup{colorlinks,linkcolor=,urlcolor=red} 

\usefonttheme{serif} % Looks like Computer Modern for non-math text -- nice!
\setbeamertemplate{navigation symbols}{} % Supress navigation symbols
\usetheme{Berlin} % Displays sections on top
\usepackage[english]{babel}
% \definecolor{links}{HTML}{2A1B81}
% \definecolor{links}{red}

\setbeamertemplate{footline}[frame number] 

\mode<presentation>
% \mode<handout>{\setbeamercolor{background canvas}{bg=black!5}}



\title{Statistical models and estimation\footnote{See last slide for copyright information.}}
\subtitle{STA431 Spring 2015}
\date{} % To suppress date


\begin{document}

\begin{frame}
  \titlepage
\end{frame}

\begin{frame}
\frametitle{Overview}
\tableofcontents
\end{frame}

\section{Models}

\begin{frame}
\frametitle{Statistical model}
\framesubtitle{Most good statistical analyses are based on a \emph{model} for the data.}

\pause

A \emph{statistical model} is a set of assertions that partly specify the probability distribution of the observable data. The specification may be direct or indirect.
 
\pause

  \begin{itemize}
    \item Let $X_1, \ldots, X_n$ be a random sample from a normal distribution with expected value $\mu$ and variance $\sigma^2$.
    
\pause

    \item For $i=1, \ldots, n$, let $Y_i = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_k x_{ik} + \epsilon_i$, where\pause
    \begin{itemize}
        \item[] $\beta_0, \ldots, \beta_k$ are unknown constants.
        \item[] $x_{ij}$ are known constants.
        \item[] $\epsilon_1, \ldots, \epsilon_n$ are independent $N(0,\sigma^2)$ random variables.
        \item[] $\sigma^2$ is an unknown constant.
        \item[] $Y_1, \ldots, Y_n$ are observable random variables.
    \end{itemize}
  \end{itemize}

 \pause

A model is not the same thing as the \emph{truth}.
\end{frame}

\begin{frame}
\frametitle{Statistical models leave something unknown}
\framesubtitle{Otherwise they are probability models} 
 \pause
\begin{itemize}
     \item The unknown part of the model for the data is called the \emph{parameter}. \pause
     \item Usually, parameters are (vectors of) numbers. \pause
     \item Usually denoted by $\theta$ or $\boldsymbol{\theta}$ or other Greek letters. \pause
     \item Parameters are unknown constants. 
\end{itemize}

\end{frame}


\begin{frame}
\frametitle{Parameter Space}

The \emph{parameter space} is the set of values that can be taken on by the parameter. 

\pause
 
  \begin{itemize}
    \item Let $X_1, \ldots, X_n$ be a random sample from a normal distribution with expected value $\mu$ and variance $\sigma^2$. 
\pause    

The parameter space is $\Theta = \{(\mu,\sigma^2): -\infty < \mu < \infty, \sigma^2 > 0\}$.

\pause

    \item For $i=1, \ldots, n$, let $Y_i = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_k x_{ik} + \epsilon_i$, where
    \begin{itemize}
        \item[] $\beta_0, \ldots, \beta_k$ are unknown constants.
        \item[] $x_{ij}$ are known constants.
        \item[] $\epsilon_1, \ldots, \epsilon_n$ are independent $N(0,\sigma^2)$ random variables.
        \item[] $\sigma^2$ is an unknown constant.
        \item[] $Y_1, \ldots, Y_n$ are observable random variables.
    \end{itemize}

\pause

The parameter space is $\Theta = \{(\beta_0, \ldots, \beta_k, \sigma^2): -\infty < \beta_j < \infty, \sigma^2 > 0\}$.
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Parameters need not be numbers}
%\framesubtitle{} 
Let $X_1, \ldots, X_n$ be a random sample from a continuous distribution with unknown distribution function $F(x)$.
\pause
\begin{itemize}
    \item The parameter is the unknown distribution function $F(x)$. \pause
    \item The parameter space is a space of distribution functions. \pause
    \item We may be interested only in a \emph{function} of the parameter, like \pause
\end{itemize}
\begin{displaymath}
    \mu = \int_{-\infty}^\infty x f(x) \, dx
\end{displaymath} \pause
\vspace{3mm}
The rest of $F(x)$ is just a nuisance parameter.
\end{frame}

\begin{frame}
\frametitle{General statement of a statistical model}
\framesubtitle{$D$ is for Data} 
{\LARGE
\begin{displaymath}
    D \sim P_\theta, ~~~ \theta \in \Theta
\end{displaymath} 
} % End size
\vspace{3mm} \pause
\begin{itemize}
    \item Both $D$ and $\theta$ could be vectors \pause
    \item For example, 
        \begin{itemize}
            \item $D = \mathbf{Y}_1, \ldots  \mathbf{Y}_n$ independent multivariate normal. \pause
            \item $\theta = (\boldsymbol{\mu,\Sigma})$. \pause
            \item $P_\theta$ is the joint distribution function of $\mathbf{Y}_1, \ldots  \mathbf{Y}_n$, with joint density \pause
        \end{itemize}
\end{itemize}
\begin{displaymath}
    f(\mathbf{y}_1, \ldots  \mathbf{y}_n) = \prod_{i=1}^n f(\mathbf{y}_i;\boldsymbol{\mu,\Sigma})
\end{displaymath} 
\end{frame}

\begin{frame}
\frametitle{Estimation}
\framesubtitle{For the model $D \sim P_\theta, ~~~ \theta \in \Theta$} 
\begin{itemize}
    \item We don't know $\theta$. \pause
    \item We never know $\theta$. \pause
    \item All we can do is guess. \pause
    \item Estimate $\theta$ (or a function of $\theta$) based on the observable data. \pause
    \item $T$ is an \emph{estimator} of $\theta$ (or a function of $\theta$): $T=T(D)$ \pause
\end{itemize}
For example, 
\begin{itemize}
    \item $D = X_1, \ldots, X_n \stackrel{i.i.d}{\sim} N(\mu,\sigma^2)$, $T = (\overline{X},S^2)$. \pause
    \item For an ordinary multiple regression model, $T=(\widehat{\boldsymbol{\beta}},MSE)$ \pause
\end{itemize}
$T$ is a \emph{statistic}, a random variable (vector) that can be computed from the data without knowing the values of any unknown parameters.
\end{frame}

\section{MOM}

\begin{frame}
\frametitle{Parameter estimation}
\framesubtitle{For the model $D \sim P_\theta, ~~~ \theta \in \Theta$} 
\begin{itemize}
    \item Estimate $\theta$ with $T=T(D)$. \pause
    \item How do we get a recipe for $T$? \pause
    \item It's good to be systematic. Lots of methods are available. \pause
    \item We will consider two: Method of moments and maximum likelihood.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Moments}
\framesubtitle{Based on a random sample like $(X_1,Y_1), \ldots, (X_n,Y_n)$}  \pause
\begin{itemize}
    \item Moments are quantities like $E\{X_i\}$, $E\{X_i^2\}$, $E\{X_iY_i\}$, $E\{W_iX_i^2Y_i^3\}$, etc. \pause
    \item \emph{Central} moments are moments of \emph{centered} random variables:  \pause
        \begin{itemize}
            \item[] $E\{(X_i-\mu_x)^2\}$ \pause
            \item[] $E\{(X_i-\mu_x)(Y_i-\mu_y)\}$ \pause
            \item[] $E\{(X_i-\mu_x)^2(Y_i-\mu_y)^3(Z_i-\mu_z)^2\}$ \pause
        \end{itemize}
    \item These are all \emph{population} moments.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Population moments and sample moments}
% \framesubtitle{Assume $X_i$ and $Y_i$ values can be observed.} 
\begin{center}
\renewcommand{\arraystretch}{1.5}
\begin{tabular}{ll} \hline
Population moment & Sample moment \\ \hline 
$E\{X_i\}$ & $\frac{1}{n}\sum_{i=1}^n X_i$  \\  \pause
$E\{X_i^2\}$ & $\frac{1}{n}\sum_{i=1}^n X_i^2$  \\ \pause
$E\{X_iY_i\}$ & $\frac{1}{n}\sum_{i=1}^n X_iY_i$  \\ \pause
$E\{(X_i-\mu_x)^2\}$ & $\frac{1}{n}\sum_{i=1}^n (X_i-\overline{X}_n)^2$  \\  \pause
$E\{(X_i-\mu_x)(Y_i-\mu_y)\}$ 
   & $\frac{1}{n}\sum_{i=1}^n (X_i-\overline{X}_n)(Y_i-\overline{Y}_n)$  \\ \pause
$E\{(X_i-\mu_x)(Y_i-\mu_y)^2\}$ 
   & $\frac{1}{n}\sum_{i=1}^n (X_i-\overline{X}_n)(Y_i-\overline{Y}_n)^2$  \\
\end{tabular}
\renewcommand{\arraystretch}{1.0}
\end{center}
\end{frame}

\begin{frame}
\frametitle{Estimation by the Method of Moments (MOM)}
\framesubtitle{For the model $D \sim P_\theta, ~~~ \theta \in \Theta$} 
\begin{itemize} \pause
    \item Population moments are a function of $\theta$. \pause
    \item Find $\theta$ as a function of the population moments. \pause
    \item Estimate $\theta$ with that function of the \emph{sample} moments. \pause
\end{itemize}
Symbolically,  \pause
\begin{itemize}
    \item Let $m$ denote a vector of population moments.  \pause
    \item $\widehat{m}$ is the corresponding vector of sample moments. \pause
    \item Find $m = g(\theta)$ \pause
    \item Solve for $\theta$, obtaining $\theta= g^{-1}(m)$. \pause
    \item Let $\widehat{\theta} = g^{-1}(\widehat{m})$. \pause
\end{itemize}
It doesn't matter if you solve first or put hats on first.
\end{frame}

\begin{frame}
\frametitle{Example: $X_1, \ldots, X_n \stackrel{i.i.d}{\sim} U(0,\theta)$}
\framesubtitle{$f(x) = \frac{1}{\theta}$ for $0<x<\theta$}  \pause
First find the moment (expected value). \pause
\begin{eqnarray*}
    E(X_i) & = & \int_0^\theta x \frac{1}{\theta} \, dx \\ \pause
           & = &  \frac{1}{\theta} \int_0^\theta x \, dx \\ \pause
           & = &  \frac{1}{\theta} \left. \frac{x^2}{2} \right|_0^\theta  \pause
             =    \frac{1}{2\theta} (\theta^2-0) \\  \pause
           & = &  \frac{\theta}{2}  \pause
\end{eqnarray*} 
\vspace{2mm}  

So $m = \frac{\theta}{2}  \pause \Leftrightarrow \theta = 2m$, \pause and $\widehat{\theta} = 2\overline{X}$. 
\end{frame}

\begin{frame}[fragile]
\frametitle{Small numerical example} \pause
%\framesubtitle{}
Let  $X_1, \ldots, X_n$ be a random sample from a uniform distribution on $(0,\theta)$. Estimate $\theta$ by the Method of Moments for the following data. Your answer is a number. Show some work.
{\footnotesize % or scriptsize
\begin{verbatim}
4.09 0.13 0.84 3.83 2.13 4.67 4.61 0.40 4.19 0.71
\end{verbatim} 
} % End size
\vspace{7mm} \pause
$\overline{X}=2.56$ so $\widehat{\theta} = 2\overline{X} = 2*2.56 = 5.12$. 
\end{frame}

\begin{frame}
\frametitle{Method of moments estimators are not unique}
\framesubtitle{What moments you use are up to you.} 
\pause
\begin{displaymath}
    E(X_i^2) = \frac{1}{\theta} \int_0^\theta x^2 \, dx = \frac{\theta^2}{3}
\end{displaymath} \pause
\vspace{2mm}
So set $m = \frac{\theta^2}{3} \Leftrightarrow \theta  = \sqrt{3m}$, and \pause
{\LARGE
\begin{displaymath}
    \widehat{\theta} = \sqrt{\frac{3}{n}\sum_{i=1}^n X_i^2}
\end{displaymath} 
} % End size
 \pause
Compared to $2\overline{X}$.
\end{frame}

\begin{frame}[fragile]
\frametitle{Compare $\widehat{\theta}_1 = 2\overline{X}$ and $\widehat{\theta}_2 = \sqrt{\frac{3}{n}\sum_{i=1}^n X_i^2}$ }
\framesubtitle{For the numerical example}  \pause
{\footnotesize % or scriptsize
\begin{verbatim}
x     4.09   0.13   0.84    3.83   2.13    4.67    4.61   0.40  4.19   0.71  
x^2  16.7281 0.0169 0.7056 14.6689 4.5369 21.8089 21.2521 0.16 17.5561 0.5041

\end{verbatim} 
} % End size

 \pause
\begin{displaymath}
    \widehat{\theta}_1 = 5.12 ~~~~~~~~ \widehat{\theta}_2 = 5.42
\end{displaymath} 
\end{frame}

\begin{frame}
\frametitle{Method of Moments estimator for normal}
\framesubtitle{Let $X_1, \ldots, X_n \stackrel{i.i.d}{\sim} N(\mu,\sigma^2)$}  \pause
From the moment-generating function or a textbook, $E(X_i)=\mu$ and $E(X_i^2) = \sigma^2 + \mu^2$. \pause  Solving for the parameters, \pause
%{\LARGE
\begin{eqnarray*}
    \mu & = & E(X_i) \\
    \sigma^2 & = & E(X_i^2) - \left( E(X_i) \right)^2 
\end{eqnarray*}  \pause
%} % End size
so
\begin{eqnarray*}
    \widehat{\mu} & = & \overline{X} \\
    \widehat{\sigma}^2 & = & \frac{1}{n}\sum_{i=1}^n X_i^2 - \overline{X}^2
\end{eqnarray*} 

\end{frame}

\begin{frame}
\frametitle{A regression example}
\framesubtitle{Independently for $i=1, \ldots, n$,}  
\begin{displaymath}
     Y_i = \beta_0 + \beta_1 X_i + \epsilon_i, \mbox{ where}
\end{displaymath} \pause
\begin{itemize}
    \item $E(X_i)=\mu_x$, $Var(X_i)=\sigma^2_x$ \pause
    \item $E(\epsilon_i)=0$, $Var(\epsilon_i)=\sigma^2_\epsilon$ \pause
    \item $X_i$ and $\epsilon_i$ are independent. \pause
    \item The distributions of $X_i$ and $\epsilon_i$ are unknown. \pause
    \item What's the parameter? \pause
    \item[]
    \item The parameter is $(\beta_0,\beta_1,F_\epsilon(\epsilon),F_x(x)\,)$. \pause
    \item We want to estimate $\beta_0$ and $\beta_1$,  \pause a \emph{function} of the parameter.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Calculate some moments}
\framesubtitle{$Y_i = \beta_0 + \beta_1 X_i + \epsilon_i$} 
 \pause
\begin{itemize}
    \item[] $E(X_i)= \mu_x$ \pause
    \item[] $Var(X_i) = \sigma^2_x$ \pause
    \item[] $E(Y_i) = \beta_0 + \beta_1 \mu_x$ \pause
    \item[] $Cov(X_i,Y_i) = \beta_1 \sigma^2_x$ \pause
\end{itemize}
Use the centering rule to get the last one: \pause
\begin{eqnarray*}
    Cov(X_i,Y_i) & = & E(\stackrel{c}{X}_i\stackrel{c}{Y}_i) \\ \pause
    & = & E\{ \stackrel{c}{X}_i (\beta_1 \stackrel{c}{X}_i + \epsilon_i) \} \\ \pause
    & = & E\{ \beta_1\stackrel{c}{X}\stackrel{2}{\vphantom{r}_i} 
          + \stackrel{c}{X}_i\epsilon_i) \} \\ \pause
    & = & \beta_1 E\{ \stackrel{c}{X}\stackrel{2}{\vphantom{r}_i} \} 
          + E\{ \stackrel{c}{X}_i\} E\{\epsilon_i\} \\ \pause
    & = &  \beta_1 \sigma^2_x
\end{eqnarray*} 
\end{frame}

\begin{frame}
\frametitle{Solve for $\beta_0$ and $\beta_1$}
\framesubtitle{Have $E(X_i)= \mu_x$, $Var(X_i) = \sigma^2_x$, 
$E(Y_i) = \beta_0 + \beta_1 \mu_x$, $Cov(X_i,Y_i) = \beta_1 \sigma^2_x$}   \pause
Putting hats on first, solve \pause
\begin{eqnarray*}
    \overline{Y} & = & \widehat{\beta}_0 + \widehat{\beta}_1 \overline{X} \\ \pause
    \widehat{\sigma}_{xy} & = & \widehat{\beta}_1 \widehat{\sigma}^2_x \pause
\end{eqnarray*}
$\Rightarrow$ \pause
\begin{eqnarray*}
    \widehat{\beta}_1 & = & \frac{\widehat{\sigma}_{xy}}{\widehat{\sigma}^2_x} \pause
    = \frac{\sum_{i=1}^n (X_i-\overline{X}_n)(Y_i-\overline{Y}_n)}
               {\sum_{i=1}^n (X_i-\overline{X}_n)^2} \mbox{ and} \\ \pause
    \widehat{\beta}_0 & = & \overline{Y} - \widehat{\beta}_1 \overline{X}
\end{eqnarray*}  \pause
These happen to be the same as the least-squares estimates.
\end{frame}

% herehere

\begin{frame}
\frametitle{Multivariate multiple regression}
\framesubtitle{Multivariate means more than one response variable} 
\begin{center}
\includegraphics[width=4in]{RegressionPath2}
\end{center} \pause
We will obtain method of moments estimation for this.
\end{frame}


\begin{frame}
\frametitle{One regression equation for each response variable}
\framesubtitle{Give the equations in scalar form.}  \pause
\begin{center}
\includegraphics[width=2in]{RegressionPath2}
\end{center} \pause

\begin{eqnarray*}
    Y_{i,1} & = & \beta_{0,1} + \beta_{1,1}X_{i,1} + \beta_{1,2}X_{i,2} + \beta_{1,3}X_{i,3}
                 + \epsilon_{i,1} \\
    Y_{i,2} & = & \beta_{0,2} + \beta_{2,1}X_{i,1} + \beta_{2,2}X_{i,2} + \beta_{2,3}X_{i,3}
                 + \epsilon_{i,2} \\
\end{eqnarray*}
\end{frame}

\begin{frame}
\frametitle{$\mathbf{Y}_i = \boldsymbol{\beta}_0 + \boldsymbol{\beta}_1 \mathbf{X}_i +   \boldsymbol{\epsilon}_i$} \pause
%\framesubtitle{} 
\begin{columns}   % Use Beamer's columns to use more of the margins!
\column{1.15\textwidth}
{\small
In scalar form, \pause
\begin{eqnarray*}
    Y_{i,1} & = & \beta_{0,1} + \beta_{1,1}X_{i,1} + \beta_{1,2}X_{i,2} + \beta_{1,3}X_{i,3}
                 + \epsilon_{i,1} \\
    Y_{i,2} & = & \beta_{0,2} + \beta_{2,1}X_{i,1} + \beta_{2,2}X_{i,2} + \beta_{2,3}X_{i,3}
                 + \epsilon_{i,2} \\
\end{eqnarray*} \pause
In matrix form, \pause
\begin{displaymath}
\begin{array}{cccccccc}
\mathbf{Y}_i  &=& \boldsymbol{\beta}_0 &+& \boldsymbol{\beta}_1 & \mathbf{X}_i &+&   \boldsymbol{\epsilon}_i  \\
& & & & & & & \\
\left( \begin{array}{c} Y_{i,1} \\  Y_{i,2}   \end{array} \right)           &=&
\left( \begin{array}{c} \beta_{0,1} \\  \beta_{0,2}   \end{array} \right)   &+&
\left( \begin{array}{ccc}
                 \beta_{1,1} & \beta_{1,2} & \beta_{1,3}  \\
                 \beta_{2,1} & \beta_{2,2} & \beta_{2,3}  \\
    \end{array} \right) &
\left( \begin{array}{c} X_{i,1} \\  X_{i,2} \\ X_{i,3} \end{array} \right)  &+&
\left( \begin{array}{c} \epsilon_{i,1} \\  \epsilon_{i,2} \end{array} \right) \\
\end{array}
\end{displaymath} 
} % End size
\end{columns}
\end{frame}


\begin{frame}
\frametitle{Statement of the model}
% \framesubtitle{These are all vectors and matrices.} 
Independently for $i=1, \ldots, n$, \pause
\begin{displaymath} 
    \mathbf{Y}_i = \boldsymbol{\beta}_0 + \boldsymbol{\beta}_1 \mathbf{X}_i +   \boldsymbol{\epsilon}_i, \mbox{ where}
\end{displaymath} \pause
{\footnotesize
\begin{itemize}
    \item $\mathbf{Y}_i$ is an $q \times 1$ random vector of observable response variables, so the regression is multivariate; there are $q$ response variables. \pause
    \item $\mathbf{X}_i$ is a $p \times 1$ observable random vector; there are $p$ explanatory variables. $E(\mathbf{X}_i) = \boldsymbol{\mu}_x$ and $V(\mathbf{X}_i) = \boldsymbol{\Phi}_{p \times p}$. The matrix $\boldsymbol{\Phi}$ is unknown. \pause
    \item $\boldsymbol{\beta}_0$ is a  $q \times 1$ matrix of unknown constants. \pause
    \item $\boldsymbol{\beta}_1$ is a  $q \times p$ matrix of unknown constants. These are the regression coefficients, with one row for each response variable and one column for each explanatory variable. \pause
    \item $\boldsymbol{\epsilon}_i$ is a $q \times 1$ random vector with expected value zero and unknown variance-covariance matrix $V(\boldsymbol{\epsilon}_i) = \boldsymbol{\Psi}_{q \times q}$.  \pause 
    \item $\boldsymbol{\epsilon}_i$ is independent of $\mathbf{X}_i$.
\end{itemize}
} % End size
\end{frame}

\begin{frame}
\frametitle{A Method of Moments estimate of $\boldsymbol{\beta}_1$}
\framesubtitle{$\mathbf{Y}_i = \boldsymbol{\beta}_0 + \boldsymbol{\beta}_1 \mathbf{X}_i +   \boldsymbol{\epsilon}_i$}  \pause
Denote the $p \times q$ matrix of (population) covariances between $\mathbf{X}_i$ and $\mathbf{Y}_i$ by  \pause
\begin{eqnarray*}
    \boldsymbol{\Sigma}_{xy} & = & C(\mathbf{X}_i,\mathbf{Y}_i) \\ \pause
    & = &  E\{\stackrel{c}{\mathbf{X}_i}\stackrel{c}{\mathbf{Y}} 
                 \stackrel{\top}{\vphantom{r}_i}     \} \\ \pause
    & = &  E\{\stackrel{c}{\mathbf{X}}_i
             \! (\boldsymbol{\beta}_1 \! \stackrel{c}{\mathbf{X}_i} 
           + \boldsymbol{\epsilon}_i)^\top \} \\ \pause          
    & = &  E\{\stackrel{c}{\mathbf{X}}_i
              (\stackrel{c}{\mathbf{X}} \stackrel{\top}{\vphantom{r}_i}
              \! \boldsymbol{\beta}_1^\top  
           + \boldsymbol{\epsilon}_i^\top) \} \\ \pause          
    & = &  E\{\stackrel{c}{\mathbf{X}}_i
              \stackrel{c}{\mathbf{X}} \stackrel{\top}{\vphantom{r}_i}
              \! \boldsymbol{\beta}_1^\top  
           + \stackrel{c}{\mathbf{X}}_i \boldsymbol{\epsilon}_i^\top \} \\ \pause  
    & = &  E\{\stackrel{c}{\mathbf{X}}_i
              \stackrel{c}{\mathbf{X}} \stackrel{\top}{\vphantom{r}_i}\}
               \boldsymbol{\beta}_1^\top   
           + E\{ \stackrel{c}{\mathbf{X}}_i \boldsymbol{\epsilon}_i^\top \} \\ \pause  
    & = &  V(\mathbf{X}_i)\boldsymbol{\beta}_1^\top + C(\mathbf{X}_i,\boldsymbol{\epsilon}_i)    \\ \pause
    & = &  \boldsymbol{\Phi} \boldsymbol{\beta}_1^\top  + \mathbf{0} \\ \pause
    & = &   \boldsymbol{\Phi} \boldsymbol{\beta}_1^\top 
\end{eqnarray*} 

\end{frame}

\begin{frame}
\frametitle{Solve for $\boldsymbol{\beta}_1$}
\framesubtitle{In terms of moments of the observable data}  \pause
%
\renewcommand{\arraystretch}{1.5}
\begin{displaymath}
\begin{array}{rrrr}
            & \boldsymbol{\Phi} \boldsymbol{\beta}_1^\top & = & \boldsymbol{\Sigma}_{xy} \\ \pause
\Rightarrow & \boldsymbol{\Phi}^{-1} \boldsymbol{\Phi} \boldsymbol{\beta}_1^\top 
            & = & \boldsymbol{\Phi}^{-1} \boldsymbol{\Sigma}_{xy} \\  \pause
\Rightarrow &  \boldsymbol{\beta}_1^\top & = & 
               \boldsymbol{\Phi}^{-1} \boldsymbol{\Sigma}_{xy} \\  \pause
\Rightarrow &  \boldsymbol{\beta}_1 
            & = &  \boldsymbol{\Sigma}_{xy}^\top (\boldsymbol{\Phi}^{-1})^\top  \\ \pause
      &     & = &  \boldsymbol{\Sigma}_{yx}\boldsymbol{\Phi}^{-1}
\end{array}
\end{displaymath}  \pause
\renewcommand{\arraystretch}{1.0}

\vspace{5mm}

%{\LARGE
Noting $\boldsymbol{\Phi} = V(\mathbf{X}_i)$ could be written $\boldsymbol{\Sigma}_x$,  \pause
have $\boldsymbol{\beta}_1 = \boldsymbol{\Sigma}_{yx}\boldsymbol{\Sigma}_x^{-1}$ 
%} % End size
\end{frame}

\begin{frame}
\frametitle{MOM estimate of $\boldsymbol{\beta}_1$ based on $\boldsymbol{\beta}_1 = \boldsymbol{\Sigma}_{yx}\boldsymbol{\Sigma}_x^{-1}$} \pause
\framesubtitle{Just put hats on.} 
{\LARGE
\begin{displaymath}
    \widehat{\boldsymbol{\beta}}_1 = \widehat{\boldsymbol{\Sigma}}_{yx} \widehat{\boldsymbol{\Sigma}}_x^{-1},
\end{displaymath}  \pause
} % End size
where
\begin{eqnarray*}
    \widehat{\boldsymbol{\Sigma}}_{yx} & = & \frac{1}{n}\sum_{i=1}^n (\mathbf{Y}_i-\overline{\mathbf{Y}}) 
                        (\mathbf{X}_i-\overline{\mathbf{X}})^\top \\ \pause
    \widehat{\boldsymbol{\Sigma}}_x & = & \frac{1}{n}\sum_{i=1}^n (\mathbf{X}_i-\overline{\mathbf{X}}) 
                        (\mathbf{X}_i-\overline{\mathbf{X}})^\top 
\end{eqnarray*}

\end{frame}

\section{MLE}

\begin{frame}
\frametitle{Maximum likelihood estimation}
\framesubtitle{A great idea from R. A. Fisher (1890-1962)}  \pause
  \begin{itemize}
    \item Given a model and a set of observed data, how should we estimate $\theta$? \pause
    \item Find the value of $\theta$ that makes the data we observed have the highest probability. \pause
    \item If the model is continuous, maximize the probability of observing data in a little region surrounding the observed data vector. \pause
    \item In either case, let $f(\mathbf{d};\theta)$ denote the joint probability density function or probability mass function evaluated at the observed data vector. \pause
    \item Maximize $L(\theta) = f(\mathbf{d};\theta)$ over all $\theta \in  \Theta$. \pause
    \item $L(\theta)$ is called the \emph{likelihood function}.
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Maximum likelihood estimation for independent random sampling}
\begin{displaymath} \pause
     D_1, \ldots, D_n  \stackrel{i.i.d.}{\sim} P_\theta, \, \theta \in \Theta.  
\end{displaymath} \pause
\begin{displaymath}
L(\theta) = \prod_{i=1}^n f(d_i;\theta), 
\end{displaymath}  \pause
where $f(d_i;\theta)$ is the density or probability mass function evaluated at $d_i$. \pause

\vspace{1mm}
\begin{itemize}
    \item Find the value of $\theta$ for which $L(\theta)$ is maximum. \pause
    \item Or equivalently, maximize $\ell(\theta) = \ln L(\theta)$. \pause
    \item The elementary approach: \pause
        \begin{itemize}
            \item Take derivatives, \pause
            \item Set derivatives to zero, \pause
            \item Solve for $\theta$,  \pause
            \item Put a hat on the answer.
        \end{itemize}
\end{itemize}
\end{frame}

\begin{frame}{Example: Coffee taste test}
A fast food chain is considering a change in the blend of coffee beans they use to make their coffee. To determine whether their customers prefer the new blend, the company plans to select a random sample of $n=100$ coffee-drinking customers and ask them to taste coffee made with the new blend and with the old blend, in cups marked ``$A$" and ``$B$." Half the time the new blend will be in cup $A$, and half the time it will be in cup $B$. Management wants to know if there is a difference in preference for the two blends.
\end{frame}

\begin{frame}{Statistical model for the taste test example} \pause
Letting $\theta$ denote the probability that a consumer will choose the new blend, treat the data $Y_1, \ldots, Y_n$ as a random sample from a Bernoulli distribution.  \pause That is, independently for $i=1, \ldots, n$,  
\begin{displaymath} \pause
    f(y_i;\theta) = \theta^{y_i} (1-\theta)^{1-y_i}
\end{displaymath}
for $y_i=0$ or $y_i=1$, and zero otherwise. 
% The conditional probability notation is not in the book (I believe).

% \vspace{10mm}

%Note that $Y = \sum_{i=1}^n Y_i$ is the number of consumers who choose the new blend. Because $Y \sim B(n,\theta)$, the whole experiment could also be treated as a single observation from a Binomial.

\end{frame}

\begin{frame}{Find the MLE of $\theta$}{Show your work}

Maximize the log likelihood. \pause
\begin{eqnarray*}
    \frac{\partial}{\partial\theta} \ln L(\theta) 
    & = & \frac{\partial}{\partial\theta} \ln\left(\prod_{i=1}^n f(y_i;\theta) \right) \\ \pause
    & = & \frac{\partial}{\partial\theta} \ln\left(\prod_{i=1}^n \theta^{y_i} (1-\theta)^{1-y_i} \right) \\ \pause
    & = & \frac{\partial}{\partial\theta} \ln\left(\theta^{\sum_{i=1}^n y_i} 
                                       (1-\theta)^{n-\sum_{i=1}^n y_i}\right) \\ \pause
    & = & \frac{\partial}{\partial\theta}\left((\sum_{i=1}^n y_i)\ln\theta + 
          (n-\sum_{i=1}^n y_i)\ln (1-\theta) \right) \\ \pause
    & = & \frac{\sum_{i=1}^n y_i}{\theta} - \frac{n-\sum_{i=1}^n y_i}{1-\theta}
\end{eqnarray*}
\end{frame}

\begin{frame}{Setting the derivative to zero,}
\begin{eqnarray*}
\frac{\sum_{i=1}^n y_i}{\theta} = \frac{n-\sum_{i=1}^n y_i}{1-\theta}   \pause
    & \Rightarrow & (1-\theta)\sum_{i=1}^n y_i = \theta (n-\sum_{i=1}^n y_i) \\ \pause
    & \Rightarrow & \sum_{i=1}^n y_i-\theta\sum_{i=1}^n y_i = n\theta - \theta\sum_{i=1}^n y_i \\ \pause
    & \Rightarrow & \sum_{i=1}^n y_i = n\theta  \\ \pause
    & \Rightarrow & \theta = \frac{\sum_{i=1}^n y_i}{n}  \pause
\end{eqnarray*}
So it looks like the MLE is the sample proportion. Carrying out the second derivative test to be sure,
\end{frame}

\begin{frame}{Second derivative test}
\begin{eqnarray*}
    \frac{\partial^2\ln \ell}{\partial\theta^2}   \pause
    & = & \frac{\partial}{\partial\theta} \left(\frac{\sum_{i=1}^n y_i}{\theta} 
                                     - \frac{n-\sum_{i=1}^n y_i}{1-\theta} \right) \\ \pause
    & = & \frac{-\sum_{i=1}^n y_i}{\theta^2} --- \frac{n-\sum_{i=1}^n y_i}{(1-\theta)^2} \\ \pause
    & = & -n\left(\frac{1-\overline{y}}{(1-\theta)^2} + \frac{\overline{y}}{\theta^2} \right)
          \pause < 0 \\
\end{eqnarray*} \pause

\vspace{5mm}

Concave down, maximum, verifying $\widehat{\theta}=\overline{y}$.
\end{frame}

\begin{frame}[fragile]{Numerical estimate} % Note use of fragile to make verbatim work.
Suppose 60 of the 100 consumers prefer the new blend. Give a point estimate the parameter $\theta$. Your answer is a number. \pause

\vspace{10mm}

\begin{verbatim}
> ybar = 60/100; ybar
[1] 0.6
\end{verbatim}
%\verb:f(x):
\end{frame}


\begin{frame}
\frametitle{Maximum likelihood for the univariate normal}
%\framesubtitle{} 
 \pause
Let $X_1, \ldots, X_n \stackrel{i.i.d}{\sim} N(\mu,\sigma^2)$. \pause

\vspace{5mm}

\begin{eqnarray*}
    \ell(\theta) & = & \ln \prod_{i=1}^n \frac{1}{\sigma\sqrt{2\pi}}
                                          e^{-\frac{1}{2}\frac{(x_i-\mu)^2}{\sigma^2}}\\ \pause
     & = & \ln\left(\sigma^{-n}(2\pi)^{-\frac{n}{2}} e^{-\frac{1}{2\sigma^2}
                                                \sum_{i=1}^n(x_i-\mu)^2}\right) \\ \pause
     & = & -n\ln\sigma - \frac{n}{2}\ln (2\pi) - \frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-\mu)^2
\end{eqnarray*} 
\end{frame}

\begin{frame}
\frametitle{Differentiate with respect to the parameters}
\framesubtitle{$\ell(\theta) = -n\ln\sigma - \frac{n}{2}\ln (2\pi) - \frac{1}{2\sigma^2}\sum_{i=1}^n(x_i-\mu)^2$}  \pause

\begin{eqnarray*}
   \frac{\partial\ell}{\partial\mu} & = & - \frac{1}{2\sigma^2}\sum_{i=1}^n2(x_i-\mu)(-1)  \pause
                                          \stackrel{set}{=} 0 \\ \pause
      & \Rightarrow & \mu= \overline{x} \\ \pause
&&\\
   \frac{\partial\ell}{\partial\sigma} & = &   \pause - \frac{n}{\sigma} 
 - \frac{1}{2} \sum_{i=1}^n(x_i-\mu)^2 (-2\sigma^{-3}) \\ \pause
& = & - \frac{n}{\sigma} + \frac{1}{\sigma^3} \sum_{i=1}^n(x_i-\mu)^2  \pause
                                         ~~ \stackrel{set}{=} ~~ 0 \\  \pause
      & \Rightarrow & \sigma^2 = \frac{1}{n}\sum_{i=1}^n(x_i-\mu)^2 
\end{eqnarray*} 
\end{frame}

\begin{frame}
\frametitle{Substituting}
%\framesubtitle{}
Setting derivaties to zero, we have obtained  \pause
\begin{displaymath}
    \mu = \overline{x} \mbox{ and } \sigma^2 = \frac{1}{n}\sum_{i=1}^n(x_i-\mu)^2, \mbox{ so}
\end{displaymath}  \pause

{\LARGE
\begin{eqnarray*}
    \widehat{\mu} & = & \overline{X} \\ \pause
    \widehat{\sigma}^2 & = & \frac{1}{n}\sum_{i=1}^n(X_i-\overline{X})^2 \\
\end{eqnarray*} 
} % End size
\end{frame}

\begin{frame}
\frametitle{Gamma Example}
%\framesubtitle{} 
Let $X_1, \ldots, X_n$ be a random sample from a Gamma distribution with parameters $\alpha>0$ and $\beta>0$  \pause
{\LARGE
\begin{eqnarray*}
     f(x;\alpha,\beta) & = & \frac{1}{\beta^\alpha \Gamma(\alpha)} 
                             e^{-x/\beta} x^{\alpha - 1} \\ \pause
&& \\
    \Theta & = & \{(\alpha,\beta): \alpha>0, \beta>0 \}  
\end{eqnarray*} 
} % End size
\end{frame}

\begin{frame}
\frametitle{Log Likelihood}
\framesubtitle{$f(x;\alpha,\beta) = \frac{1}{\beta^\alpha \Gamma(\alpha)} 
                             e^{-x/\beta} x^{\alpha - 1}$} 
\begin{eqnarray*}
\ell(\alpha,\beta) 
    &=& \ln \prod_{i=1}^n \frac{1}{\beta^\alpha \Gamma(\alpha)} 
                 e^{-x_i/\beta} x_i^{\alpha - 1} \nonumber \\ \pause
    &=& \ln \left( \beta^{-n\alpha} \, \Gamma(\alpha)^{-n}
                 \exp(-\frac{1}{\beta}\sum_{i=1}^n x_i) 
                 \left(\prod_{i=1}^n x_i \right)^{\alpha-1} 
            \right)  \\ \pause
    &=& -n\alpha\ln\beta -n\ln\Gamma(\alpha) - \frac{1}{\beta}\sum_{i=1}^n x_i
        + (\alpha - 1) \sum_{i=1}^n \ln x_i    
\end{eqnarray*} 
\end{frame}

\begin{frame}
\frametitle{Differentiate with respect to the parameters}
\framesubtitle{$\ell(\theta) = -n\alpha\ln\beta -n\ln\Gamma(\alpha) - \frac{1}{\beta}\sum_{i=1}^n x_i + (\alpha - 1) \sum_{i=1}^n \ln x_i$}  \pause

\begin{eqnarray*}
   \frac{\partial\ell}{\partial\beta} & \stackrel{set}{=} & 0~~ 
                 \pause   \Rightarrow ~~ \alpha\beta = \overline{x}  \\ \pause
&&\\
   \frac{\partial\ell}{\partial\alpha} & = &   -n\ln\beta 
            -n  \frac{\partial}{\partial\alpha} \ln \Gamma(\alpha) + \sum_{i=1}^n \ln x_i   \\ \pause
& = & \sum_{i=1}^n \ln x_i    -n\ln\beta   - n\frac{\Gamma^\prime(\alpha)}{\Gamma(\alpha)}
         \pause  ~~\stackrel{set}{=}~~0        \\ 
\end{eqnarray*} 
\end{frame}

\begin{frame}
\frametitle{Solve for $\alpha$} 
%\framesubtitle{} 
{\LARGE
\begin{displaymath}
    \sum_{i=1}^n \ln x_i    -n\ln\beta   - n\frac{\Gamma^\prime(\alpha)}{\Gamma(\alpha)} = 0
\end{displaymath}  \pause
} % End size

\vspace{8mm}

where 
\begin{displaymath}
    \Gamma(\alpha) = \int_0^\infty e^{-t}t^{\alpha-1} \, dt.
\end{displaymath} \pause

\vspace{8mm}

Nobody can do it.
\end{frame}

\begin{frame}
\frametitle{Maximize the likelihood numerically with software}
\framesubtitle{Usually this is in high dimension} 
\begin{center}
\includegraphics[width=4in]{Likelihood}
\end{center} \pause

\begin{itemize}
    \item It's like trying to find the top of a mountain by walking uphill blindfolded. \pause
    \item You might stop at a local maximum. \pause
    \item The starting place is very important. \pause
    \item The final answer is a number (or vector of numbers). \pause
    \item There is no explicit formula for the MLE. 
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{There is a lot of useful theory}
\framesubtitle{Even without an explicit formula for the MLE}  
\begin{center}
\includegraphics[width=3in]{Likelihood}
\end{center} \pause

\begin{itemize}
    \item MLE is asymptotically normal. \pause
    \item Variance of the MLE is deeply related to the curvature of the log likelihood at the MLE. \pause
    \item The more curvature, the smaller the variance. \pause
    \item The variance of the MLE can be estimated from the curvature (using the Fisher Information). \pause
    \item Basis of tests and confidence intervals. 
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Comparing MOM and MLE} \pause
%\framesubtitle{} 
  \begin{itemize}
    \item Sometimes they are identical, sometimes not. \pause
    \item If the model is right they are usually close for large samples. \pause
    \item Both are asymptotically normal. \pause
    \item Estimates of the variance are well known for both. \pause
    \item Small variance of an estimator is good. \pause
    \item As $n \rightarrow \infty$, nothing can beat the MLE. \pause
    \item Except that the MLE depends on a very specific distribution. \pause
    \item And sometimes the dependence matters. 
  \end{itemize}
\end{frame}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{frame}
\frametitle{Copyright Information}

This slide show was prepared by  \href{http://www.utstat.toronto.edu/~brunner}{Jerry Brunner},
Department of Statistics, University of Toronto. It is licensed under a 
\href{http://creativecommons.org/licenses/by-sa/3.0/deed.en_US}
     {Creative Commons Attribution - ShareAlike 3.0 Unported License}. Use any part of it as you like and share the result freely. The \LaTeX~source code is available from the course website:
\href{http://www.utstat.toronto.edu/~brunner/oldclass/431s15} {\footnotesize \texttt{http://www.utstat.toronto.edu/$^\sim$brunner/oldclass/431s15}}

\end{frame}


\end{document}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{}
%\framesubtitle{} 
  \begin{itemize}
    \item 
    \item 
    \item 
  \end{itemize}
\end{frame}

% \stackrel{c}{\mathbf{X}} \stackrel{\top}{\vphantom{r}_i}