% These are my slides for my ASA talk about my research.

\documentclass{slides}

%\documentclass{seminar}
%\documentclass[article]{seminar}
%\special{landscape}
%\twoup[1]
\pagestyle{empty}
%\onlyslides{28,29}

\input epsf.tex

\usepackage{amsbsy}

\newcommand{\bp}{\mbox{\boldmath $\phi$}}
\newcommand{\bT}{\mbox{\boldmath $\Theta$}}

\begin{document}

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

\begin{slide}
\begin{center}
{\bf Bayesian Deconvolution of Seismic Array Data} 
\end{center}
\begin{center}
{\bf for Ripple-Fired Explosions}
\end{center}
\begin{center}
Eric A. Suess \\
\end{center}
\begin{center}
Advisor: Robert Shumway \\
\end{center}
\end{slide}

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

\begin{slide}
{\bf Outline:}

\begin{enumerate}
%\item
%History
\item
The Problem
\item
Model
\item
Bayesian Statistics
\item
Bayesian Deconvolution of Seismic Array Data
\item
Results
\begin{itemize}
\item
Simulated Data
\item
Real Data
\end{itemize}
\item
Conclusions and Future Work
\end{enumerate}

\end{slide}

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

\begin{slide}
{\bf Treaties:}

\begin{enumerate}
\item
1963 Limited Nuclear Test Ban Treaty (LTBT)
\item
1968 Non-Proliferation of Nuclear Weapons Treaty (NPT)
\item
1974 Threshold Test Ban Treaty (TTBT)
\item
1976 Peaceful Nuclear Explosions Treaty (PNET)
\item
1996 Comprehensive Test Ban Treaty (CTBT)
\end{enumerate}

\end{slide}

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

\begin{slide}
{\bf Background:}

Much of the focus in the past has been on distinguishing possible
\underline{nuclear explosions} from \underline{earthquakes}.  

Currently, since the testing treaties have put limitations on the permissible sizes of the nuclear
explosions, other smaller seismic events such as \underline{industrial
mining explosions} have become of interest in the discrimination problem.

The work presented here is related to distinguishing
low-level nuclear explosions from ripple-fired mining explosions that
are on the same seismic level.
\end{slide}

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

\begin{slide}
{\bf The Problem:}
\begin{itemize}
\item
Monitoring seismic events at Regional Distances for low-level
nuclear tests.
\item
Other seismic sources need to be ruled out.
\begin{itemize}
\item
Ripple-Fired Mining Explosions
\end{itemize}
\end{itemize}
\end{slide}

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

\begin{slide}
{\bf Ripple-Fired Explosions:}

This is a mining technique in which explosions of single devices
(or groups of devices) are \mbox{detonated} in succession.
\end{slide}

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

\begin{slide}
{\bf Monitoring:}

Arrays of receivers are put in place at \mbox{Regional} Distances and
seismic data is continually \mbox{collected}.  Seismic disturbances that
are above the baseline noise of the area are investigated.
\end{slide}

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

\begin{slide}
{\bf Model:}
\begin{eqnarray}
y_k(t) = s_k(t)+\sum_{j=1}^m a_js_k(t-j) + \varepsilon_k(t) \nonumber
\end{eqnarray}

{\it Amplitudes} are distributed according to a random
Bernoulli-Gaussian model.
{\small (ref. Cheng, Chen and Li 1996) }
\begin{eqnarray}
p(a_j|\eta) \sim (1-\eta)I(a_j=0)+\eta TN(\mu_\alpha,\sigma^2_\alpha)I(a_j > 0)
\nonumber
\end{eqnarray}

{\it Signal and path effects} follow an $AR(3)$ model. 
{\small (ref. Dargahi-Noubary 1995, Tj{\o}stheim 1975)}
\nonumber
\begin{eqnarray}
s_k(t) + \phi_1s_k(t-1) + ...+ \phi_ps_k(t-p) = e_k(t)
\nonumber
\end{eqnarray}

$e_k(t)$ i.i.d $N(0,\sigma^2)$ and define the precision
$\tau=\frac{1}{\sigma^2}$.

$\varepsilon_k(t)$ i.i.d. $N(0,c\sigma^2)$, $ c =1/SNR$, $c \ge 0$ is
fixed.

\end{slide}

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

\begin{slide}
Truncated Normal:

\begin{eqnarray}
p(a_j) = c^{-1}(\mu_{\alpha},\sigma_{\alpha}^2)\frac{1}{\sqrt{2\pi}\sigma_{\alpha}}
\exp\left\{-\frac{(a_j-\mu_{\alpha})^2}{\sigma_{\alpha}^2}\right\}I(a_j>0) \nonumber
\end{eqnarray}

%\noindent
where
\begin{eqnarray}
c^{-1}(\mu_{\alpha},\sigma_{\alpha}^2) = \int_0^\infty \frac{1}{\sqrt{2\pi}\sigma_{\alpha}}
\exp\left\{-\frac{(x-\mu_{\alpha})^2}{\sigma_{\alpha}^2}\right\}
\nonumber
\end{eqnarray}
\end{slide}

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

\begin{slide}
{\bf Bayesian Statistics:}

Model:
\begin{eqnarray}
p({\bf Y}|\bT) \nonumber
\end{eqnarray}

Prior distribution:
\begin{eqnarray}
p(\bT) \nonumber
\end{eqnarray}

Joint distribution:
\begin{eqnarray}
p({\bf Y},\bT)=p({\bf Y}|\bT)p(\bT) \nonumber
\end{eqnarray}

Posterior distribution:
\begin{eqnarray}
p(\bT|{\bf Y})&=&\frac{p({\bf Y}|\bT)p(\bT)}{\int p({\bf Y}|\bT)
p(\bT)d\bT} \nonumber \\
&& \nonumber\\
&\propto&p({\bf Y}|\bT)p(\bT) \nonumber
\end{eqnarray}

\end{slide}

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

\begin{slide}
{\bf Gibbs Sampler:} {\small(ref. Gelfand and Smith 1990)}

$\bT = (\theta_1,\theta_2)$

$p(\bT|{\bf Y}) = p((\theta_1,\theta_2)|{\bf Y})$

%Want the marginals,
%\begin{eqnarray}
%p(\theta_1|{\bf Y}) \nonumber
%\end{eqnarray}
%\begin{eqnarray}
%p(\theta_2|{\bf Y}) \nonumber
%\end{eqnarray}

%Can get the conditional marginals,
%\begin{eqnarray}
%p(\theta_1|{\bf Y}, \theta_2)\nonumber
%\end{eqnarray}
%\begin{eqnarray}
%p(\theta_2|{\bf Y},\theta_1)\nonumber
%\end{eqnarray}

%\end{slide}

%\begin{slide}

Given $\theta_1^{(0)}$ and $\theta_2^{(0)}$ for $h = 1,...,Reps$

\begin{enumerate}
\item
Sample $\theta_1^{(h)}$ from $p(\theta_1|{\bf Y}, \theta_2^{(h-1)})$
\item
Sample $\theta_2^{(h)}$ from $p(\theta_2|{\bf Y}, \theta_1^{(h)})$
\item
Set $h = h+1$ and go to 1.
\end{enumerate}

\begin{eqnarray}
(\theta_1^{(BurnIn+1)},\theta_2^{(BurnIn+1)}),...,(\theta_1^{(Reps)},\theta_2^{(Reps)})\nonumber
\end{eqnarray}
%a ``random sample'' from the target distribution
%$p((\theta_1,\theta_2)|Y)$.

%\begin{eqnarray}
%\theta_1^{(Burn)},...,\theta_1^{(Reps)} \nonumber
%\end{eqnarray}
%a ``random sample'' from the
%marginal $p(\theta_1|Y)$.

\end{slide}

\begin{slide}

$\bT^{(1)},...,\bT^{(Reps)}$ are realizations of a stationary Markov
Chain, with transition probability from $\bT^{(h-1)}$ to $\bT^{(h)}$,
\begin{eqnarray}
T(\bT^{(h-1)},\bT^{(h)}) =
p(\theta_1|{\bf Y},\theta_2^{(h-1)})p(\theta_2|{\bf Y},\theta_1^{(h)})
\nonumber
\end{eqnarray}

By Ergodic theory, we can calculate \mbox{estimates} of say
$\theta_1$ by
\begin{eqnarray}
\frac{1}{Reps}\sum_h \theta_1^{(i)} \mathop{\longrightarrow}\limits^{a.s.}
E[\theta_1|{\bf Y}], \nonumber \\
\nonumber \\
Reps \mathop{\rightarrow} \infty. \nonumber
\end{eqnarray}

\end{slide}

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

\begin{slide}

{\bf Priors:}

%$\eta$ = the {\it probability} of seeing $a_j \ne 0$ at time $j$.

\begin{eqnarray}
\eta \sim BETA(\beta_1,\beta_2)
\nonumber
\end{eqnarray}

%The $AR$ coefficients $\bp = [ \phi_1, \phi_2,...,\phi_p]^\prime$.

\begin{eqnarray}
\bp \sim N_p(\bp_0,\Sigma_0)
\nonumber
\end{eqnarray}


%\nonumber
%\end{eqnarray}

%$\bp_0$ and $\Sigma_0$ are predetermined.

%$\tau$ = the {\it precision} of the noise $e_k(t)$.

\begin{eqnarray}
\tau \sim GAMMA(\gamma_1,\gamma_2)
\nonumber
\end{eqnarray}

%$\gamma_1 > 0$ and $\gamma_1 > 0$ are predetermined.
\end{slide}

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

\begin{slide}
{\bf The Model:}

\begin{eqnarray}
y_k(t) = s_k(t)+\sum_{j=1}^m a_js_k(t-j) + \varepsilon_k(t) \nonumber
\end{eqnarray}

{\bf The parameter set:}

$\bT = \{\eta,{\bf a},\bp,\tau,{\bf S}\}$

{\bf Hyperparamters:}

$\beta_1,\beta_2,\mu_\alpha,\sigma^2_\alpha,\bp_0,\Sigma_0,
\gamma_1,\gamma_2$

{\bf Fixed parameters:}

c, m

\end{slide}

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

\begin{slide}
{\bf Likelihood:}

${\bf Y} = [{\bf y}_1,...,{\bf y}_q]^\prime$

${\bf y}_k = [y_k(1),...,y_k(n)]^\prime$

\begin{eqnarray}
p({\bf Y} | \bT)&=&\prod_{k=1}^q p({\bf y}_k | \bT)
\nonumber \\
&=& \prod_{k=1}^q (2\pi)^{-n/2}
\left(\frac{\tau}{c}\right)^{n/2} \exp \left\{-\frac{\tau}{2c}
\sum_t\varepsilon^2_k(t)\right\}  \nonumber
\end{eqnarray}

\end{slide}

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

\begin{slide}
{\bf Overall Prior:}

\begin{eqnarray}
p(\bT)&=&p(\eta,{\bf a}, \bp,\tau,{\bf S}) \nonumber \\
\nonumber \\
&=&p(\eta)p({\bf a}|\eta)p(\bp)p(\tau)p({\bf S}|\bp,\tau)
\nonumber
\end{eqnarray}

{\bf Joint Density:}

\begin{eqnarray}
p({\bf Y},\bT)=p({\bf Y}|\bT)p(\bT) \nonumber
\end{eqnarray}

{\bf Joint Posterior:}

\begin{eqnarray}
p(\bT|{\bf Y}) \propto p({\bf Y}|\bT)p(\bT)
\nonumber
\end{eqnarray}

\end{slide}

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

\begin{slide}
{\bf Conditional Marginal Posterior Distributions:}

\begin{eqnarray}
p(\eta|{\bf Y},rest)\sim beta(\beta_1^*,\beta_2^*) \nonumber
\end{eqnarray}

For fixed $j = 1,...,m$
{\small
\begin{eqnarray}
p(a_j|{\bf Y},rest)\sim(1-\eta_j) I(a_j=0) +
\eta_jTN(\mu_{a_j},\sigma_{a_j}^{2})I(a_j > 0) \nonumber
\end{eqnarray}  }
\begin{eqnarray}
p(\bp|{\bf Y},rest)\sim N_p(\bp_*,\Sigma_*) \nonumber
\end{eqnarray}
\begin{eqnarray}
p(\tau|{\bf Y},rest) \sim gamma(\gamma_1^*,\gamma_2^*) \nonumber
\end{eqnarray}

For fixed $i=1,...,n$ and $k = 1,...,q$
\begin{eqnarray}
p(s_k(i)|{\bf Y},rest)\sim N(\mu_{s_k(i)},\sigma^2_{s_k(i)})
\nonumber
\end{eqnarray}

\end{slide}

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

\begin{slide}

\begin{eqnarray}
\eta|{\bf Y},rest \sim beta(\beta_1^*,\beta_2^*)
\nonumber
\end{eqnarray}

\begin{eqnarray}
\beta_1^* = m - n_a + \beta_1,
\nonumber
\end{eqnarray}

\begin{eqnarray}
\beta_2^* = n_a + \beta_2
\nonumber
\end{eqnarray}

\end{slide}

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

\begin{slide}

$a_j|{\bf Y},rest \sim$ Bernoulli-Gaussian 

\begin{eqnarray*}
\mu_{a_j} = \sigma_{a_j}^2\left[\frac{\mu_{\alpha}}{\sigma_{\alpha}^2} +
\left(\frac{\tau}{c}\right)\sum_k\sum_t\varepsilon_k^*(t[-j])s_k(t-j)\right]
\end{eqnarray*}

\begin{eqnarray*}
\sigma_{a_j}^{-2}=\frac{1}{\sigma_{\alpha}}
+\left(\frac{\tau}{c}\right)\sum_k\sum_ts^2_k(t-j)
\end{eqnarray*}

\begin{eqnarray*}
\eta_j = \frac{\eta}{\eta +
(1-\eta)\left[\frac{c(\mu_{a_j},\sigma_{a_j}^2)}{c(\mu_{\alpha},\sigma_{\alpha}^2)}\right] \left(\frac{\sigma_{a_j}}{\sigma_{\alpha}}\right)\exp\left\{\frac{1}{2}
\left[\frac{\mu^2_{a_j}}{\sigma^2_{a_j}} -
\frac{\mu^2_\alpha}{\sigma^2_{\alpha}}\right]\right\} }
\end{eqnarray*}

\end{slide}

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

\begin{slide}

\begin{eqnarray}
\bp|{\bf Y},rest \sim N_p(\bp_*,\Sigma_*) 
\nonumber
\end{eqnarray}

\begin{eqnarray*}
\bp_* = \Sigma_*\left[-\tau\sum_k\sum_t s_k(t)\tilde{\bf s}_k(t) + \Sigma_0^{-1}\bp_0\right]
\end{eqnarray*}

\begin{eqnarray*}
\Sigma_*^{-1} = \tau \sum_k\sum_t \tilde{\bf s}_k(t)\tilde{\bf s}^\prime_k(t) +\Sigma_0^{-1}
\end{eqnarray*}

where
\begin{eqnarray*}
\tilde{\bf s}_k(t) = [s_k(t-1),...,s_k(t-p)]^\prime
\end{eqnarray*}

\end{slide}

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

\begin{slide}

\begin{eqnarray}
\tau|{\bf Y},rest \sim gamma(\gamma_1^*,\gamma_2^*)
\nonumber
\end{eqnarray}

\begin{eqnarray*}
\gamma_1^* = q(2n - l - p)/2 + \gamma_1
\end{eqnarray*}

\begin{eqnarray*}
\gamma_2^* = \frac{1}{2c}\sum_k\sum_t\varepsilon_k^2(t) + \frac{1}{2}\sum_k\sum_t e_k^2(t)
\end{eqnarray*}

\end{slide}

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

\begin{slide}
$s_k(t)|{\bf Y},rest \sim$ Normal 

\begin{eqnarray*}
\mu_{s_k(i)} = \sigma^2_{s_k(i)}\left[\left(\frac{\tau}{c}\right)
\sum_ta_{t-i}\varepsilon^\prime(t[-i])
-\tau\sum_t\phi_{t-i}e_k^\prime(t[-i])\right]
\end{eqnarray*}

\begin{eqnarray*}
\sigma_{s_k(i)}^{-2} = \left(\frac{\tau}{c}\right)\sum_t a^2_{t-i}
+\tau\sum_t \phi^2_{t-i}
\end{eqnarray*}

\end{slide}

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

\begin{slide}
{\bf Steps To Perform The Gibbs Sampler:}

Given the initial values:
\begin{eqnarray}
\left\{ \eta^{(0)},{\bf a}^{(0)},\bp^{(0)},\tau^{(0)},{\bf S}^{(0)}
\right\} \nonumber
\end{eqnarray}
for $h = 1$ to Reps:

\begin{enumerate}

\item
Sample $\eta^{(h)}$ from a {\it beta}.

\item
Sample $a_j^{(h)}$, $j=1,...,m$, from a Bernoulli-Gaussian.

\item
Sample $\bp^{(h)}$ from a p-variate {\it normal}.

\item
Sample $\tau^{(h)}$ from a {\it gamma}.

\item
Sample $s_k(i)^{(h)}$, $i={p+1},...,n$, $k=1,...,q$, from a {\it normal}.

\item
Set $h = h+1$ and go to Step 1.
\end{enumerate}

\end{slide}

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

\begin{slide}
{\bf Conclusion and Further Work:}

\begin{enumerate}
\item
one
\item
two
\end{enumerate}

\end{slide}

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

\begin{slide}
{\bf References:}

Cheng, Q., Chen, R., and Li,T.-H. (1996) Simultaneous Wavelet Estimation and Deconvolution of Reflection Seismic Signals \emph{IEEE Transactions On Geoscience and Remote Sensing} {\bf 34}, 377-384.

Dargahi-Noubary. (1992) Stochastic Modeling and Identification of Seismic Records
Based on Established Deterministic Formulations. \emph{Journal of Time
Series Analysis }
{\bf 16}, 201-220.

Tj{\o}stheim, D. (1975) Autoregressive Representation of Seismic P-wave Signals with
an Application to the Problem of Short-Period Discriminants \emph{Geophys. J. R. astr. Soc.}
{\bf 43}, 269-291.

\end{slide}

\begin{slide}

Gelfand, A.E. and Smith, A.F.M. (1990) Sampling-Based Approaches to Calculating Marginal Densities \emph{Journal of the American Statistical Association}
{\bf 85}, 398-409.

Pearson, D.C., Strump, B.W., and Anderson, D.P. (1996) Physical Constraints on Mining Explosions
\emph{http://www.geology.smu.edu/~dpa-www/papers/srl/SRL.html}

\end{slide}

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

\end{document}