%% ================================================================
\section{Problem statement}

In 1986, the space shuttle Challenger exploded during takeoff, killing the
seven astronauts aboard.  It was determined that the explosion was the result
of an O-ring failure.  The flight accident was believed to be caused by the
unusually cold weather ($31^\circ$F) at the time of the launch.  The past
O-ring failure data along with temperature at launch time are given below (in
increasing order of temperature) for 18 prior flights.  The flight numbers
denote the (unimportant) time order of launch.  The numbers 0 and 1 indicate no
O-ring faiilure and O-ring failure, respectively.

%14  9 23 10  1  5 13 15  4  3  8 17  2
% 1  1  1  1  0  0  0  0  0  0  0  0  1
%53 57 58 63 66 67 67 67 68 69 70 70 70
%
%11  6  7 16 21 19 22 12 20 18
% 1  0  0  0  1  0  0  0  0  0
%70 72 73 75 75 76 76 78 79 81

%14 1 53
% 9 1 57
%23 1 58
%10 1 63
% 1 0 66
% 5 0 67
%13 0 67
%15 0 67
% 4 0 68
% 3 0 69
% 8 0 70
%17 0 70
% 2 1 70
%11 1 70
% 6 0 72
% 7 0 73
%16 0 75
%21 1 75
%19 0 76
%22 0 76
%12 0 78
%20 0 79
%18 0 81

%14  9 23 10  1  5 13 15  4  3  8 17  2  11  6  7 16 21 19 22 12 20 18
% 1  1  1  1  0  0  0  0  0  0  0  0  1   1  0  0  0  1  0  0  0  0  0
%53 57 58 63 66 67 67 67 68 69 70 70 70  70 72 73 75 75 76 76 78 79 81

$$
\begin{tabular}{|r|r|r|}
\hline
Flight \# & Failure $f_i$ & Temp. $t_i$ in ${}^\circ$F \\
\hline
\hline 14 & 1 & 53 \\
\hline  9 & 1 & 57 \\
\hline 23 & 1 & 58 \\
\hline 10 & 1 & 63 \\
\hline  1 & 0 & 66 \\
\hline  5 & 0 & 67 \\
\hline 13 & 0 & 67 \\
\hline 15 & 0 & 67 \\
\hline  4 & 0 & 68 \\
\hline  3 & 0 & 69 \\
\hline  8 & 0 & 70 \\
\hline 17 & 0 & 70 \\
\hline  2 & 1 & 70 \\
\hline 11 & 1 & 70 \\
\hline  6 & 0 & 72 \\
\hline  7 & 0 & 73 \\
\hline 16 & 0 & 75 \\
\hline 21 & 1 & 75 \\
\hline 19 & 0 & 76 \\
\hline 22 & 0 & 76 \\
\hline 12 & 0 & 78 \\
\hline 20 & 0 & 79 \\
\hline 18 & 0 & 81 \\
\hline
\end{tabular}
$$

Estimate the probability of O-ring failure at temperature $31^\circ$F.

%% ----------------------------------------------------------------
\subsection{Suggested model}

Let $F$ denote the failure status (response
variable) and $T$ the temperature in ${}^\circ$F at launch time (explanatory
variable).  Use the \emph{logistic regression model}:
\begin{align*}
	P(F=1 \mid T=t; \alpha, \beta) &= \frac{e^{\alpha+\beta t}}
		{1+e^{\alpha+\beta t}} := p_{\alpha,\beta}(t) \\
	P(F=0 \mid T=t; \alpha, \beta) &= \frac{1}
		{1+e^{\alpha+\beta t}} = 1 - p_{\alpha,\beta}(t).
\end{align*}
which may be summarized as
\begin{align}
\label{eqn:P_alpha_beta}
	P(F=f \mid T=t; \alpha, \beta)
	&= \frac{1-f + e^{\alpha+\beta t}f}{1+e^{\alpha+\beta t}} \\
	&= \frac{1+(e^{\alpha+\beta t}-1)f}{1+e^{\alpha+\beta t}}.
\end{align}
Note that one may express the model as
$$
	\log\left(\frac{p_{\alpha,\beta}(t)}{1-p_{\alpha,\beta}(t)}\right)
	= \alpha + \beta t,
$$
hence the name \emph{logistic regression}.

For the $n=23$ independent observations $\vecf=(f_1, \ldots, f_{23})$ and
$\vect=(t_1, \ldots, t_{23})$, the likelihood function of the observations,
given $(\alpha, \beta)$ is
\begin{align}
\label{eqn:ell_alpha_beta}
	\ell(\vecf \mid \vect ; \alpha, \beta) =
	\prod_{i=1}^{23}
		p_{\alpha,\beta}(t_i)^{f_i}
		(1-p_{\alpha,\beta}(t_i))^{1-f_i}.
\end{align}

%% ----------------------------------------------------------------
\subsection{Suggested method}
\label{sec:suggested_method}

Use a suitable prior distribution of $(\alpha, \beta)$.  For example, you may
let $\eta := \ln(-\beta)$ and assume $\alpha$ and $\eta$ to be independent
random variables with $\alpha$ having the normal distribution $N(\alphahat,
\sigma_1^2)$ and  $\eta$   having the normal distribution $N(\etahat,
\sigma_2^2)$.  Observe that in the present context $\beta < 0$.  (Note:  If you
choose this prior, I will supply the values of the MLEs $\alphahat$ and
$\betahat$ from Project 1, as well as the bootstrap variances $\sigma_1^2$ and
$\sigma_2^2$, if you like.)

(Later note:  Take $\alphahat = 15$ and $\betahat=-0.23$, so $\etahat=-1.47$.
Use $\sigma_1^2 =6$  and $\sigma_2^2 =0.12$.)

\begin{itemize}

\item[(a)] Compute the posterior distribution $\pi(\alpha,\eta \mid \vecf,
\vect)$ of $\alpha$ and $\eta$.

\item [(b)] Using Markov chain Monte Carlo, generate a sequence of $(\alpha_k,
\eta_k)$ --- say, compute 100,000 of them and take every 50th (discarding the
0th and 50th).  This will be approximately an i.i.d. sample from
$\pi(\alpha,\eta \mid \vecf, \vect)$.  For each $(\alpha_k, \eta_k)$ pair,
compute $p_{\alpha_k,\eta_k} (31)$ and $p_{\alpha_k,\eta_k} (t)$ for some
other temperature(s), e.g. 65.

\item[(c)] Form a histogram of the $p_{\alpha_k,\eta_k}(31)$'s.  Locate the
10th percentile $\delta_{0.1}$ of the histogram for the posterior distribution
in (i) for the failure probability at $t=31^\circ$F.  Observe that, in the
Bayesian paradigm, the statistician has a 90\% confidence that the probability
of failure at $31^\circ$F is at least $\delta_{0.1}$.  Repeat for desired $t$
values other than $t=31$.

\end{itemize}

%% ================================================================
\section{Posterior distribution}

%% ----------------------------------------------------------------
\subsection{Background}
\label{sec:bayesian_background}

Recall that given a random variable $X$ with parameter $\theta$,
the Bayesian paradigm consists of the following.

\begin{itemize}

\item Start with a \emphidx{prior distribution} $p(\theta)$, which is
the (perhaps ad hoc) knowledge of the parameter before the observations.

\item Write down the \emphidx{sampling distribution} $\phi(x\mid \theta)$.

\item Take $\vecx = (x_1, \ldots, x_n)$ i.i.d. samples from the distribution
	with unknown parameter $\theta$.  Write down the joint sampling
	distribution $\Phi(\vecx \mid \theta)$.

\item The \emphidx{posterior distribution} of the parameter $\theta$ after the
observations have been taken is
	\begin{align}
	\label{eqn:bayes_theta_x}
		\pi(\theta \mid \vecx) = \frac
			{\Phi(\vecx \mid \theta) \; p(\theta)}
			{m(\vecx)}
	\end{align}
	where the \emphidx{marginal distribution} is
	$$
		m(\vecx) = \int \Phi(\vecx \mid \theta) \pi(\theta) \, d\theta.
	$$
	Mnemonically, equation \ref{eqn:bayes_theta_x} becomes
	$$
		\textrm{posterior} =
		\frac{ \textrm{sampling} \times \textrm{prior} }
		{ \textrm{marginal} }.
	$$

\end{itemize}

%% ----------------------------------------------------------------
\subsection{Re-expression in terms of $(\alpha,\eta)$}

The first task is to convert the logistic model from $(\alpha,\beta)$
parameters to $(\alpha,\eta)$ parameters.  Since $\eta = \ln(-\beta)$, we have
$\beta=-e^\eta$ and $e^{\beta t} = e^{-e^\eta t}$.  Then the probability for
failure at a given temperature (equation \ref{eqn:P_alpha_beta}) becomes
\begin{align}
\label{eqn:P_alpha_eta}
	P(F=f \mid T=t; \alpha, \eta)
	&= \frac{1+(e^{\alpha+\beta t}-1)f}{1+e^{\alpha+\beta t}} \\
	&= \frac{1+(e^{\alpha-e^\eta t}-1)f}{1+e^{\alpha-e^\eta t}}.
\end{align}
Using abusive notation, we continue to write
$$
	p_{\alpha,\eta}(t) = P(F=1 \mid T=t; \alpha, \eta)
$$
and
\begin{align}
\label{eqn:ell_alpha_eta}
	\ell(\vecf \mid \vect ; \alpha, \eta) =
	\prod_{i=1}^{23}
		p_{\alpha,\eta}(t_i)^{f_i}
		(1-p_{\alpha,\eta}(t_i))^{1-f_i}.
\end{align}

%% ----------------------------------------------------------------
\subsection{Posterior distribution of $(\alpha,\eta)$}

Following the recipe of section \ref{sec:bayesian_background}, with parameters
$(\alpha, \eta)$ and observations $(\vecf, \vect)$, we have the following.

\begin{itemize}

\item The prior distributions of $\alpha$ and $\eta$ are independent
	$\sim N(\alphahat, \sigma_1^2)$ and $N(\etahat, \sigma_2^2)$.  Thus, the
	joint prior distribution of $\alpha$ and $\eta$ factors as $q_1(\alpha) \,
	q_2(\eta)$.  Explicitly, we have
	$$
		q_1(\alpha) = \frac{1}{\sqrt{2\pi \sigma_1^2}}
			e^{-\frac{1}{2}\left(\frac{\alpha-\alphahat}{\sigma_1}\right)^2}
	$$
	and
	$$
		q_2(\eta) = \frac{1}{\sqrt{2\pi \sigma_2^2}}
			e^{-\frac{1}{2}\left(\frac{\eta-\etahat}{\sigma_2}\right)^2}.
	$$

\item We do not quite have the joint sampling distribution
$\Phi(\vecf, \vect \mid \alpha, \eta)$.
Rather, we have the likelihood function
$\ell(\vecf \mid \vect, \alpha, \eta)$.  Let's solve for the former.
Recall from elementary probability that
$$
	\ell(\vecf \mid \vect, \alpha, \eta)
	=
	\frac{\Phi(\vecf, \vect \mid \alpha, \eta)}
		{\int \Phi(\vecf, \vect \mid \alpha, \eta) \, d\vecf}.
$$
Since each of the $n$ $f_i$'s takes on the values 0,1 we have
\begin{align*}
	\ell(\vecf \mid \vect, \alpha, \eta)
	&=
	\frac{\Phi(\vecf, \vect \mid \alpha, \eta)}
		{\sum_{f_1=0,1} \ldots
		\sum_{f_n=0,1}
		\Phi(\vecf, \vect \mid \alpha, \eta)} \\
	&=
	\frac{\Phi(\vecf, \vect \mid \alpha, \eta)}
		{\psi(\vect \mid \alpha, \eta)}
\end{align*}
where $\psi$ is the (unknown) distribution of $t$ given $\alpha$ and $\eta$.
Then the sampling distribution is
\begin{align*}
	\Phi(\vecf, \vect \mid \alpha, \eta)
	&=
	\ell(\vecf \mid \vect, \alpha, \eta)
	\; \psi(\vect \mid \alpha, \eta).
\end{align*}

\item The posterior distribution of the parameters $\alpha$ and $\eta$ after
the observations $\vecf=(f_1, \ldots, f_n)$ and $\vect=(t_1, \ldots, t_n)$ have
been taken is
	\begin{align*}
		\pi(\alpha, \eta \mid \vecf, \vect) = \frac
			{\Phi(\vecf, \vect \mid \alpha, \eta) \; 
				q_1(\alpha) \; q_2(\eta)}
			{m(\vecf, \vect)}
	\end{align*}
	where the marginal distribution $m(\vecf, \vect)$ is the integral of the
	numerator with respect to $\alpha$ and $\eta$.
	Writing this in terms of known expressions, we obtain
	\begin{align}
	\label{eqn:bayes_alpha_eta}
		\pi(\alpha, \eta \mid \vecf, \vect) = \frac{
				\ell(\vecf \mid \vect, \alpha, \eta) \;
				\psi(\vect \mid \alpha, \eta) \;
				q_1(\alpha) \;
				q_2(\eta)}
			{m(\vecf, \vect)}.
	\end{align}

\end{itemize}

%% ================================================================
\section{MCMC sampling}

As suggested in section \ref{sec:suggested_method}, we will generate a sequence
of $(\alpha_k, \eta_k)$ pairs sampled from the posterior distribution which is
given by equation \ref{eqn:bayes_alpha_eta}.  This is a bivariate MCMC,
as follows:

\begin{itemize}

\item Start with $\alpha_1 = \alphahat$ and $\eta_1 = \etahat$.

\item Given $\alpha_k$ and $\eta_k$, propose
	$$
	\alpha_{k+1}' \sim N(\alpha_k, \sigma_1^2)
	$$
	and
	$$
	\eta_{k+1}' \sim N(\eta_k, \sigma_2^2).
	$$

\item Compute
	$$
		r_1 = \min\left\{
			\frac
				{\pi_1(\alpha_{k+1}' \mid \eta_k, \vecf, \vect)}
				{\pi_1(\alpha_k      \mid \eta_k, \vecf, \vect)},
			1
		\right\}
	$$
	and
	$$
		r_2 = \min\left\{
			\frac
				{\pi_2(\eta_{k+1}' \mid \alpha_k, \vecf, \vect)}
				{\pi_2(\eta_k      \mid \alpha_k, \vecf, \vect)},
			1
		\right\}.
	$$

\item With probability $r_1$, accept the change
$$\alpha_{k+1} = \alpha_{k+1}';$$
with probability $1-r_1$, reject the change by setting
$$\alpha_{k+1} = \alpha_k.$$
Likewise, accept/reject $\eta_{k+1} = \eta_{k+1}'$ or $\eta_k$ with
probabilities $r_2$ and $1-r_2$, respectively.

\item For each $(\alpha_k, \eta_k)$ pair (appropriately downsampled as
described below), compute the realization of the random variable $p_{\alpha_k,
\eta_k}(31)$.  Sort these $p$'s and take the 10th percentile value.

\end{itemize}

We compute
$$
	\pi_1(\alpha \mid \eta, \vecf, \vect)
	=
	\frac{\pi(\alpha, \eta \mid \vecf, \vect)}
		{\int \pi(\alpha, \eta \mid \vecf, \vect) \; d\alpha}
$$
and
$$
	\pi_2(\eta \mid \alpha, \vecf, \vect)
	=
	\frac{\pi(\alpha, \eta \mid \vecf, \vect)}
		{\int \pi(\alpha, \eta \mid \vecf, \vect) \; d\eta}.
$$
Note that we only need the ratios of $\pi_1$'s and $\pi_2$'s in the
accept/reject stage of the MCMC algorithm.  Thus
\begin{align*}
	r_1
	&= \min\left\{
		\frac
			{\pi_1(\alpha_{k+1}' \mid \eta_k, \vecf, \vect)}
			{\pi_1(\alpha_k      \mid \eta_k, \vecf, \vect)},
		1
	\right\} \\
%
	&= \min\left\{
		\frac
			{\pi(\alpha_{k+1}', \eta_k \mid \vecf, \vect)}
			{\pi(\alpha_k,      \eta_k \mid \vecf, \vect)}
		\frac
			{\int \pi(\alpha, \eta_k \mid \vecf, \vect) \; d\alpha}
			{\int \pi(\alpha, \eta_k \mid \vecf, \vect) \; d\alpha},
		1
	\right\} \\
%
	&= \min\left\{
		\frac
			{\pi(\alpha_{k+1}', \eta_k \mid \vecf, \vect)}
			{\pi(\alpha_k,      \eta_k \mid \vecf, \vect)},
		1
	\right\}.
\end{align*}
Expanding the numerator and denominator as in equation
\ref{eqn:bayes_alpha_eta}, we have (with the marginal $m$ canceling
immediately)
\begin{align*}
	r_1
%
	&= \min\left\{
		\frac{
			\ell(\vecf \mid \vect, \alpha_{k+1}', \eta_k) \;
		} {
			\ell(\vecf \mid \vect, \alpha_k, \eta_k) \;
		}
		\frac{
			\psi(\vect \mid \alpha_{k+1}', \eta_k) \;
		}{
			\psi(\vect \mid \alpha_k, \eta_k) \;
		}
		\frac{
			q_1(\alpha_{k+1}')
		}{
			q_1(\alpha_k)
		}
		\frac{
			q_2(\eta_k)
		}{
			q_2(\eta_k)
		},
		1
		\right\} \\
%
	&= \min\left\{
		\frac{
			\ell(\vecf \mid \vect, \alpha_{k+1}', \eta_k) \;
		} {
			\ell(\vecf \mid \vect, \alpha_k, \eta_k) \;
		}
		\frac{
			\psi(\vect \mid \alpha_{k+1}', \eta_k) \;
		}{
			\psi(\vect \mid \alpha_k, \eta_k) \;
		}
		\frac{
			q_1(\alpha_{k+1}')
		}{
			q_1(\alpha_k)
		},
		1
		\right\}.
\end{align*}
Assuming (reasonably!) that the distribution of temperatures is independent
of the failure-modeling parameters $\alpha$ and $\eta$, this simplifies
to
\begin{align*}
	r_1
	&= \min\left\{
		\frac{
			\ell(\vecf \mid \vect, \alpha_{k+1}', \eta_k) \;
		} {
			\ell(\vecf \mid \vect, \alpha_k, \eta_k) \;
		}
		\frac{
			q_1(\alpha_{k+1}')
		}{
			q_1(\alpha_k)
		},
		1
		\right\}.
\end{align*}
Last, the scaling constant in $q_1$ cancels so we are left with
\begin{align*}
	r_1
	&= \min\left\{
		\frac{
			\ell(\vecf \mid \vect, \alpha_{k+1}', \eta_k) \;
		} {
			\ell(\vecf \mid \vect, \alpha_k, \eta_k) \;
		}
		\frac{
			e^{-\frac{1}{2}\left(\frac{\alpha_{k+1}'-\alphahat}{\sigma_1}
				\right)^2}
		}{
			e^{-\frac{1}{2}\left(\frac{\alpha_k-\alphahat}{\sigma_1}\right)^2}
		},
		1
		\right\} \\
%
	&= \min\left\{
		\frac{
			\ell(\vecf \mid \vect, \alpha_{k+1}', \eta_k) \;
		} {
			\ell(\vecf \mid \vect, \alpha_k, \eta_k) \;
		}
		\;
		e^{-\frac{1}{2\sigma_1^2}\left(
			(\alpha_{k+1}'- \alphahat)^2 -
			(\alpha_k     - \alphahat)^2
			\right)
		},
		1
		\right\}.
\end{align*}

Likewise,
\begin{align*}
	r_2
	&= \min\left\{
		\frac{
			\ell(\vecf \mid \vect, \alpha_k, \eta_{k+1}') \;
		} {
			\ell(\vecf \mid \vect, \alpha_k, \eta_k) \;
		}
		\;
		e^{-\frac{1}{2\sigma_2^2}\left(
			(\eta_{k+1}'- \etahat)^2 -
			(\eta_k     - \etahat)^2
			\right)
		},
		1
		\right\}.
\end{align*}

It is well known (see for example B.~Berg's \emph{Markov Chain Monte Carlo
Simulations And Their Statistical Analysis}) that the first few MCMC iterates
should be discarded, while the implicit Markov chain approaches its stationary
distribution:  this stage is known as \emph{thermalization}.  As well, one may
choose to retain only every $M$th (e.g. 50th) $(\alpha_k, \eta_k)$ pair, since
adjacent samples tend to be correlated.

%% ================================================================
\section{Results}

%% ----------------------------------------------------------------
\subsection{Raw MCMC data}

For informative purposes only, here we plot the sequences
of $\alpha_k$, $\eta_k$, and $p_{\alpha_k, \eta_k}(t)$ generated during
MCMC runs.

\begin{center}
\includegraphics{31.eps}

\includegraphics{65.eps}

\includegraphics{80.eps}
\end{center}

The variation in the $p_{\alpha_k,\eta_k}(t)$ values is lost at these scales.
Here are plots of $k$ vs.  $p_{\alpha_k,\eta_k}(t)$.

\begin{center}
\includegraphics{31p.eps}

\includegraphics{65p.eps}

\includegraphics{80p.eps}
\end{center}

%% ----------------------------------------------------------------
\subsection{Failure probability}

Here we collect the 10th percentile of $p_{\alpha_k, \eta_k}(t)$
for various values of $t$, in tabular and graphical form.

$$
\begin{tabular}{|c|c|}
\hline Temperature $t$ in ${}^\circ$F & $p(t)$ \\
\hline
\hline  31 & 0.9982 \\
\hline  35 & 0.9958 \\
\hline  40 & 0.9892 \\
\hline  45 & 0.9698 \\
\hline  50 & 0.9173 \\
\hline  55 & 0.7924 \\
\hline  60 & 0.5635 \\
\hline  65 & 0.2901 \\
\hline  70 & 0.1161 \\
\hline  75 & 0.0368 \\
\hline  80 & 0.0118 \\
\hline  85 & 0.0032 \\
\hline  90 & 0.0008 \\
\hline  95 & 0.0002 \\
\hline 100 & 0.0001 \\
\hline
\end{tabular}
$$

% 31 35 40 45 50 55 60 65 70 75 80 85 90 95 100

\begin{center}
\includegraphics{failprob.eps}
\end{center}

%% ================================================================
\section{Conclusion}

Based on the logistic regression model and Bayesian parameter estimation used
here, it is clear that O-ring failure was almost certain at a launch
temperature of $31^\circ$F.
