\documentclass[11pt]{article}
\usepackage{latexsym}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{epsfig}


\newcommand{\handout}[5]{
  \noindent
  \begin{center}
  \framebox{
    \vbox{
      \hbox to 5.78in { {#1} \hfill #2 }
      \vspace{4mm}
      \hbox to 5.78in { {\Large \hfill #5  \hfill} }
      \vspace{2mm}
      \hbox to 5.78in { {\em #3 \hfill #4} }
    }
  }
  \end{center}
  \vspace*{4mm}
}

\newcommand{\lecture}[5]{\handout{#1}{#2}{#3}{#4}{#5}}

\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{observation}[theorem]{Observation}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{claim}[theorem]{Claim}
\newtheorem{fact}[theorem]{Fact}
\newtheorem{assumption}[theorem]{Assumption}

% 1-inch margins, from fullpage.sty by H.Partl, Version 2, Dec. 15, 1988.
\topmargin 0pt
\advance \topmargin by -\headheight
\advance \topmargin by -\headsep
\textheight 8.9in
\oddsidemargin 0pt
\evensidemargin \oddsidemargin
\marginparwidth 0.5in
\textwidth 6.5in

\parindent 0in
\parskip 1.5ex

%\renewcommand{\baselinestretch}{1.25}

\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Document header: change the lecture number, date, scribe name, and
% lecture topic.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\lecture{Statistical Techniques in Robotics (16-831, F10)}
{Lecture\#18 (Thursday October 21)}
{Lecturer: Drew Bagnell}
{Scribe: Heather Justice \footnotemark}
{Bayesian Linear Regression Part 2 \& Gaussian Properties}

% Be sure to acknowledge all previous scribes.
\footnotetext{Some content adapted from previous scribes: Laura Lindzey, Hans Pirnay, Siddharth Mehrotra}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Document body: Scribe your notes here.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Bayesian Linear Regression}

\subsection{Problem Formulation}

Recall from the previous lecture our formulation of the Bayesian Linear Regression problem, with the corresponding graphical model shown in Figure~\ref{fig-BLRModel}. The variables $x_0,...,x_t$ are our input feature vectors. The weight vector is $\theta \sim \mathcal{N}(\mu_\theta, \Sigma_\theta)$. Our prediction is based on:
\begin{eqnarray*}
y_t \sim \mathcal{N}(\theta^Tx_t, \sigma^2) = \theta^Tx_t + \epsilon \text{, where } \epsilon \sim \mathcal{N}(0, \sigma^2)
\end{eqnarray*}

We need to compute $P(\theta|D)$, where $D$ is all data $x_0,...,x_t$ and $y_0,...,y_t$, so that we can make a prediction for $y_{t+1}$ by conditioning on $\theta$ as follows:
\begin{eqnarray*}
P(y_{t+1}|D) = \int_{\theta} P(y_{t+1}|x_{t+1},\theta) P(\theta|D) \delta\theta
\end{eqnarray*}

\begin{figure}[h]
\centering
   \includegraphics[scale=.5]{BLRModel.png}
\caption{Bayesian Linear Regression Model}
\label{fig-BLRModel}
\end{figure}

\subsection{Update Step}

Recall from the previous lecture that we can write the weight vector in the natural parameterization:
\begin{eqnarray*}
\theta \sim \tilde{\mathcal{N}}(J_\theta, P_\theta) = 
\frac{1}{z}\exp{\left( -\frac{1}{2}\theta^T P \theta + J^T \theta \right)}
\end{eqnarray*}

We do the update step as follows (as derived in the previous lecture):
\begin{eqnarray*}
J_{t+1} &=& J_t + \frac{1}{\sigma^2} y_t x_t \\
P_{t+1} &=& P_t + \frac{1}{\sigma^2} x_t x_t^T
\end{eqnarray*}

Notice that, due to the properties of Gaussians, the precision is always increasing for each new data point (in other words, the variance never increases). Also notice that the update for $P$ has no dependence on $y_t$ (the variance depends only on inputs and not outputs); this model does not account for surprises in the data! Gaussians don't really have a good model for outliers. In this model, if an outlier is observed, that will essentially shift the mean, and the algorithm will become more certain about that new mean. In practice, you may want to have some sort of wrapper function to detect and discard outliers from the data first (imagine if the data contained an outlier $10\sigma$ from the mean!).

\subsection{Transfer to Moment Parameterization}

Given the transfer rules to the moment parameterization
\begin{eqnarray*}
\Sigma &=& P^{-1}\\
\mu &=& P^{-1}J
\end{eqnarray*}
the moment parameterization after $N$ timesteps is then
\begin{eqnarray*}
\Sigma_{\theta | D} &=& \left[ \Sigma_0^{-1} + \sum_{i=1}^N \frac{x_ix_i^T}{\sigma^2} \right]^{-1}\\
\mu_\theta &=& \Sigma_{\theta | D} \cdot \left[ \sum_{t=1}^N \frac{y_tx_t}{\sigma^2} \right] \\
&=& \left[ \Sigma_0^{-1} + \sum_{i=1}^N \frac{x_ix_i^T}{\sigma^2} \right]^{-1} \cdot \left[ \sum_{t=1}^N \frac{y_tx_t}{\sigma^2} \right]
\end{eqnarray*}

Watch out for the different uses of the sigma symbol (variance versus summation)! Also note that this assumes the initial $\mu_0 = 0$; if $\mu_0$ is nonzero, then the $\mu_\theta$ update will be more complicated.

The complexity of this this computation is essentially cubed in the number of features. More precisely, $\mathcal{O}(F^3+TF^2)$, where $F$ is the number of features and $T$ is the number of data points.

\subsection{Making Predictions}

Given all data $D$ up to timestep $t$ and $x_{t+1}$, the probability of an observation $\tilde{y}_{t+1}$ is
\begin{eqnarray}
p(\tilde{y}_{t+1}|x_{t+1},D) &=& \int p(\tilde{y}_{t+1}|x_{t+1},D,\theta)\cdot p(\theta|D)d\theta\\
&=& \int p(\tilde{y}_{t+1}|x_{t+1},\theta)\cdot p(\theta,D)d\theta
\end{eqnarray}

We might expect that $E[y_{t+1}|D] = \mu_\theta^Tx_{t+1}$ since we know that $E[\theta]=\mu_\theta$ and we want $\theta$ such that $y_{t+1}=\theta^Tx_{t+1}$. More formally, we can compute:
\begin{eqnarray*}
E_{P(\theta|D)}[y_{t+1}|x_{t+1}] &=& E[\theta^Tx_{t+1}+\epsilon] \\
&=& E[\theta^Tx_{t+1}] + E[\epsilon] \\
&=& E[\theta]^Tx_{t+1} + 0 \text{ (since that $\epsilon$ is independent of $\theta$)} \\
&=& \mu_\theta^Tx_{t+1}
\end{eqnarray*}

Notice that inference is exact for this model!

\section{Properties of Gaussian Parameterizations}

There are two common parameterizations for Gaussians, the moment parameterization and the natural parameterization. It is often most practical to switch back and forth between representations, depending on which calculations are needed. The moment parameterization is more convenient for visualization (simply draw a Gaussian centered around the mean with width determined by the variance), calculating expected value, and computing marginals. The natural parameterization is more convenient for multiplying Gaussians and for conditioning on known variables.

\subsection{Moment Parameterization}

Recall that the moment parameterization of a Gaussian is:
\begin{eqnarray}
\mathcal{N}(\mu,\Sigma) = p(\theta) &=& \frac{1}{z}\exp\left( -\frac{1}{2}\left(\theta -\mu\right)^T\Sigma^{-1}\left(\theta -\mu\right)\right)
\end{eqnarray}

Given:
\begin{eqnarray*}
\mathcal{N}\left(\left[\begin{array}{c}\mu_1\\\mu_2\end{array}\right],\left[\begin{array}{cc}\Sigma_{11}&\Sigma_{12}\\\Sigma_{21}&\Sigma_{22}\end{array}\right]\right)
\end{eqnarray*}
\textbf{Marginal:} computing $p(x_2)$
\begin{eqnarray*}
\mu_2^{\mathrm{marg}} &=& \mu_2\\
\Sigma_2^{\mathrm{marg}} &=& \Sigma_{22}
\end{eqnarray*}
\textbf{Conditional:} computing $p(x_1|x_2)$
\begin{eqnarray*}
\mu_{1|2} &=& \mu_1 + \Sigma_{12}\Sigma_{22}^{-1}\left(x_2 -\mu_2\right)\\
\Sigma_{1|2} &=& \Sigma_{11} - \Sigma_{12}\Sigma_{22}^{-1} \Sigma_{21}
\end{eqnarray*}

\subsection{Natural Parameterization}

Recall that the natural parameterization of a Gaussian is:
\begin{eqnarray}
\tilde{\mathcal{N}}(J,P) = \tilde{p}(\theta) &=& \frac{1}{z}\exp\left( J^T\theta-\frac{1}{2}\theta^{T}P\theta\right)
\end{eqnarray}

Given:
\begin{eqnarray*}
\mathcal{N}\left(\left[\begin{array}{c}J_1\\J_2\end{array}\right],\left[\begin{array}{cc}P_{11}&P_{12}\\P_{21}&P_{22}\end{array}\right]\right)
\end{eqnarray*}
\textbf{Marginal:} computing $p(x_2)$
\begin{eqnarray*}
J_2^{\mathrm{marg}} &=& J_2 - P_{21}P_{11}^{-1}J_1\\
P_2^{\mathrm{marg}} &=& P_{22} - P_{21}P_{11}^{-1}P_{12}
\end{eqnarray*}
\textbf{Conditional:} computing $p(x_1|x_2)$
\begin{eqnarray*}
J_{1|2} &=& J_1 - P_{12}x_2\\
P_{1|2} &=& P_{11}
\end{eqnarray*}

\section{Introduction to the Gauss-Markov Model}

Consider $X_1,X_2,...,X_t,X_{t+1}$ to be the state variables and $Y_1,Y_2,...,Y_t,Y_{t+1}$ to be the sequence of corresponding observations. As in Hidden Markov models, conditional independencies (see Figure~\ref{fig-GaussMarkovModel}) dictate that past and future states are decorrelated given the current state, $X_t$ at time $t$. For example, if we know what $X_2$ is, then no information
about $X_1$ can possibly help us to reason about what $X_3$ should be. 
\begin{figure}[h]
\centering
   \includegraphics[width=12cm, angle=0, trim=0.25in 2.5in 0.25in 2.5in]{GaussMarkovModelDiagram.pdf}
\caption{The Independence Diagram of a Gauss-Markov Model}
\label{fig-GaussMarkovModel}
\end{figure}

Note that we want to be able to preserve Gaussian properties in this model, so our motion and observation models must be linear systems.

\textbf{Prior:} \begin{eqnarray*}x_1 \sim \mathcal{N}(\mu_0, \Sigma_0)\end{eqnarray*}

\textbf{Motion Model:} \begin{eqnarray*}x_{t+1} = Ax_t+\epsilon \text{, where } \epsilon \sim \mathcal{N}(0, Q)\end{eqnarray*}
(Note that $A$ could be indexed by time, but for simplicity we will not consider that case here.)

\textbf{Observation Model:} \begin{eqnarray*}y_{t+1} = Cx_{t+1} + \delta \text{, where } \delta \sim \mathcal{N}(0, R)\end{eqnarray*}

Ultimately we want to be able to compute $P(x_t | y_1,...,y_t)$. For example, we might want to compute the robot's location ($x_t$) at some time $t$ given the robot's observations $(y_1,...,y_t)$.

\subsection*{Introduction to the Lazy Gauss-Markov Filter}

At time step $t$, given that we have already computed $P(x_t|D) \sim \mathcal{N}(\mu_t, \Sigma_t)$, the prediction (or motion update) is easiest to do in the moment form of the Gaussian.

The form of this update is $P(x_{t+1}|D) \sim \mathcal{N}(A\mu_t, ???)$. Note that we will be deriving the variance ``???'' in a future homework. We can easily prove the expected value as follows:
\[E[X_{t+1}] = E[AX_t + \epsilon]\]
\[\Rightarrow E[X_{t+1}] = E[AX_t] + E[\epsilon]\]
since variance of $\epsilon$ is 0,
\[\Rightarrow E[X_{t+1}] = AE[X_t] = A\mu_t\] 

Further details about the Lazy Gauss-Markov Filter will be explored in the next lecture!

\end{document}
