\chapter{Methodology}
\label{ch:method}


\section{Preliminaries}

In this assignment we seek coefficients for the model

\begin{equation}
y = \sum^n_{i=1}x_i b_i + \epsilon
\end{equation}

where $n<m$, $x_i$'s are the predictors, $y$ is the response variable, and $\epsilon$ is the measurement noise.  In this problem, given $k$ independent measurements, denoted by $\mathbf{y}$, $\mathbf{X}$, and $\epsilon$.  To solve this problem, the will formulate a Bayesian solution to the joint estimation of $\{n,b_1, \ldots, b_n\}$.

\section{Definitions}

To setup a Bayesian formulation, define a joint posterior density of the type:

\begin{equation*}
f(n,\mathbf{b}_n|\mathbf{y}) \propto f(\mathbf{y},n|\mathbf{b}_n) f(\mathbf{b}_n|n)f(n)
\end{equation*}

Let the notation $X_n = X(:,1:n)$ and $\mathbf{b}_n = \{b_1, \ldots, b_n\}$.  In this project, will use the following terms

\begin{description}
\item The \underline{likelihood function} is given by
\begin{equation}
f(\mathbf{y},n|\mathbf{b}_n) = \left(\frac{1}{\sqrt{2\pi \sigma_0^2}} \right)^k e^{-\dfrac{1}{2\sigma^2_0}||\mathbf{y} - \mathbf{X_n}\mathbf{b_n} ||^2}
\end{equation}

\item The \underline{prior} on $\mathbf{b}_n$ given $n$ is:
\begin{equation}
f(\mathbf{b}_n|n) = \left(\frac{1}{\sqrt{2\pi \sigma_p^2}} \right)^k e^{-\dfrac{1}{2\sigma^2_p}||\mathbf{b}_n - \mu_b ||^2}
\end{equation}

\item The \underline{prior} on n is uniform: 
\begin{equation}
f(n) = \dfrac{1}{m}
\end{equation}

\end{description}

\section{Sampling from the Posterior}

The method of the project will involve using the Reversible Jump Monte Carlo sampling (RJMCMC) algorithm for sampling from the posterior.  
\vspace{.25cm}
\\
Let $(n,b_n)$ be the current samples from the posterior.  The following is the algorithm used to apply RJMCMC:
\begin{description}
\item[(a)] Select a candidate number, $n*$ from the probability $f(n)$.

\item[(b)] If $n* \leq n$, generate a random vector $u \sim N(0, \sigma_rI_{n*})$.  The candidate coefficient vector is given by:

\begin{equation*}
b_{n*} = 
\begin{bmatrix}
\mathbf{b}_n \\
0
\end{bmatrix}
+
\begin{bmatrix}
\mathbf{u}_1 \\
\mathbf{u}_2
\end{bmatrix}, \hspace{.5cm}
\mathbf{u}= 
\begin{bmatrix}
\mathbf{u}_1 \\
\mathbf{u}_2
\end{bmatrix}
\end{equation*}


Compute the likelihoods:
\begin{eqnarray}
h_1(\mathbf{u}) &=& \left(\frac{1}{\sqrt{2\pi \sigma_r^2}} \right)^{n*} e^{-\dfrac{1}{2\sigma^2_r}||\mathbf{u} ||^2} \nonumber
\\
h_2(\mathbf{u_1}) &=& \left(\frac{1}{\sqrt{2\pi \sigma_r^2}} \right)^{n} e^{-\dfrac{1}{2\sigma^2_r}||\mathbf{u_1} ||^2} \nonumber
\end{eqnarray}

\item[(c)] If $n* < n$, generate a random vector $u_1 \sim N(0, \sigma_rI_{n*})$.  The candidate coefficient vector is given by: 

\begin{equation*}
b_{n*} = 
\mathbf{b}_n^1
+
\mathbf{u}_1,  \hspace{.5cm}
\mathbf{b}_n = 
\begin{bmatrix}
\mathbf{b}^1_n \\
\mathbf{b}^2_n
\end{bmatrix}
\end{equation*}

and form 
$\mathbf{u} = 
\begin{bmatrix}
\mathbf{u}_1 \\
\mathbf{b}^2_n
\end{bmatrix}
$
\vspace{.25cm}
\\
Compute the likelihoods:
\begin{eqnarray}
h_2(\mathbf{u}) &=& \left(\frac{1}{\sqrt{2\pi \sigma_r^2}} \right)^{n} e^{-\dfrac{1}{2\sigma^2_r}||\mathbf{u} ||^2} \nonumber
\\
h_1(\mathbf{u_1}) &=& \left(\frac{1}{\sqrt{2\pi \sigma_r^2}} \right)^{n*} e^{-\dfrac{1}{2\sigma^2_r}||\mathbf{u_1} ||^2} \nonumber
\end{eqnarray}

\item[(d)]  Compute the acceptance-rejection function:

\begin{eqnarray}
\rho &=& \min\left(1,\dfrac{f(n^{*},\mathbf{b}_{n^{*}}|\mathbf{y})h_2}{f(n,\mathbf{b}_n|\mathbf{y})h_1}\right) = \min \left(1,\dfrac{f(\mathbf{y}|n^{*},\mathbf{b}_{n^{*}})f(\mathbf{b}_{n^{*}}|n^{*})h_2}{f(\mathbf{y}|n^{*},\mathbf{b}_n)f(\mathbf{b}_n|n)h_1}\right) \nonumber
\\
\; \nonumber
\\
&=& \min\left(1,\dfrac{e^{(E_1-E_2)}(2\pi\sigma^2_p)^{\frac{n-n^{*}}{2}}e^{-\frac{1}{2\sigma^2_p}(||\mathbf{b_{n^{*}}-\mu_b}||^2-||b_n - \mu_b||^2)}h_2}{h_1}\right) \nonumber
\end{eqnarray}

Where $E_1 = \frac{1}{2\sigma^2_0}||\mathbf{y}-\mathbf{X}_{n^{*}}\mathbf{b}_{n^{*}}||^2$ and 
$E_2 =\frac{1}{2\sigma^2_0}||\mathbf{y}-\mathbf{X}_{n}\mathbf{b}_{n}||^2$


\item[(e)]  If $U \sim U[0,1]$ is less than $\rho$ then set $(n,\mathbf{b}_n)= (n^{*},\mathbf{b}_{n^{*}})$.  Else, return to Step (a).
\end{description}
