\documentclass[11pt,twocolumn]{article}
\usepackage{times,amsmath,amsfonts}
\usepackage[dvips]{graphicx,graphics}



\setlength{\textheight}{10in} \setlength{\textwidth}{7.1in}
\setlength{\columnsep}{0.3in} \setlength{\topmargin}{0in}
\setlength{\headheight}{-0.3in} \setlength{\headsep}{0in}
\setlength{\oddsidemargin}{-0.6in}  % Centers text.


% Add the period after section numbers.  Adjust spacing.
%\newcommand{\Section}[1]{\vspace{-5pt}\section{\hskip -1em.~~#1}\vspace{-3pt}}
%\newcommand{\SubSection}[1]{\vspace{-2pt}\subsection{\hskip -1em.~~#1}
%        \vspace{-3pt}}
\newcommand{\bqn}{\begin{eqnarray*}}
\newcommand{\eqn}{\end{eqnarray*}}
\newcommand{\bq}{\begin{eqnarray}}
\newcommand{\eq}{\end{eqnarray}}

\begin{document}

\title{Stat 992: Lecture 18\\
Multiple comparison}
\author{Moo K. Chung \tt{mchung@stat.wisc.edu}}
\maketitle \thispagestyle{empty}

%\begin{figure}
%\centering
%\renewcommand{\baselinestretch}{1}
%\includegraphics[scale=0.8]{homework05-1.eps}
%\end{figure}

{\bf Problem 6.} Is the variance of signal smaller after
smoothing?\\

{\bf Theorem.} Assume the covariance function of field $Y$ to be
isotopic, i.e. $R_Y(x,y) = f(\|x-y\|)$ and $f$ decreasing. Then
$${\bf Var} [K_{\sigma}*Y(x)] \leq {\bf Var} Y(x) \mbox{ for each } x \in \mathbb{R}^n.$$
{\em Proof.} (Shijie Tang) Note that $R_Y(x',y') = f(\|x'-y'\|)
\leq f(0) = R_Y(x',x')$. Let us denote ${\bf Var} Y(x) = f(0) =
\varsigma^2$. Isotropic covariance function basically implies
stationary uniform variance field.


The covariance function $R$ of $K_{\sigma}*Y(x)$ is given by
 \bqn R(x,y)&=&\int K_{\sigma}(x-x')K_{\sigma}(y-y')R_Y(x',y')\;dx'dy'\\
&\leq& \int K_{\sigma}(x-x')K_{\sigma}(y-y')\varsigma^2\;dx'dy'\\
&\leq& \varsigma^2.\eqn By letting $x=y$, we prove the statement.
$\Box$

While investigating the above theorem, I realized that we did not
actually prove the iterated kernel identity for truncated and
normalized kernel! Note that
$$\widetilde{K}_{\sigma}*\widetilde{K}_{\sigma} \neq\widetilde{K}_{\sqrt{2}\sigma}.$$
But our numerical implementation of kernel smoothing is based on
iterated kernel smoothing so how we grantee the convergence to the
integral version of kernel smoothing?

 {\bf Theorem.} $$K_{\sigma}*K_{\sigma}(x) \leq
\widetilde{K}_{\sigma}*\widetilde{K}_{\sigma}(x) \leq
\frac{1}{\alpha^2(\Omega)} K_{\sigma}*K_{\sigma}(x), x \in
\Omega$$
 where
$\alpha(\Omega) = \int_{\Omega} K_{\sigma}(x) \; dx.$ Obviously
this can be generalized into $k$-iterated truncated normalized
kernel smoothing so that
$$K_{\sigma}^{(k)}(x) \leq
\widetilde{K}_{\sigma}^{(k)}(x) \leq \frac{1}{\alpha^k(\Omega)}
K_{\sigma}^{(k)}(x), x \in \Omega$$



For the proper choice of $\sigma$ and $\Omega$, we can make
$\alpha = 0.9999$. The usual number of iterations is less than
100-200 in practical application. That gives $\alpha^{100}
=0.9900$ and $\alpha^{200} =0.9802$. For $\alpha=0.999$,
$\alpha^{100} =0.9048$ and $\alpha^{200} = 0.8186$ Hence,
decreasing the side of bandwidth and increasing
the number of iterations would perform better.\\
 {\em Proof.}
\bqn K_{\sigma}(x) = \alpha \widetilde{K}_{\sigma}(x) +
K_{\sigma}(x)I_{\Omega^c}(x).\eqn Apply convolution to the above
equation. \bqn K_{\sigma}*K_{\sigma}(x) &\geq& \alpha K_{\sigma}
*\widetilde{K}_{\sigma}(x)\\
&\geq&
\alpha^2\widetilde{K}_{\sigma}(x)*\widetilde{K}_{\sigma}(x). \eqn
Trivially $ K_{\sigma}(x) \leq \widetilde{K}_{\sigma}(x), x \in
\Omega$ since $0< \alpha < 1$. $\Box$ This proposition shows that
the iterated Gaussian kernel smoothing can approximate the
integral version of kernel smoothing with less than one percent
accuracy if one desires in practical applications.

\begin{enumerate}

 \item {\em Multiple comparison.} We have an image $Y$ of size
$95 \times 68$ where each pixel intensity is assumed to follow
Gaussian. Suppose we have model
$$Y(x) =\mu(x) + \epsilon(x)$$
where $\epsilon$ is a zero mean unit variance Gaussian field. In
imaging, one of the most important problem is that of detection,
which can be stated as the problem of identifying the region of
statistically significant region. So it can be formulated as an
inference problem $$H_0: \mu(x)=0 \mbox{ for all } x \mbox{ vs. }
H_1:\mu(x)
> 0 \mbox { for some } x.$$
Note that the null hypothesis is based on multiple hypotheses. Let
$H_0^j:\mu(x_j) = 0, x_j \in \Omega$ and assume $\Omega =
\{x_1,\cdots, x_m\}$. Then
$$H_0 = \bigcap_{j=1}^m H_0^j.$$
We also have corresponding alternate hypothesis $H_1^j: \mu(x_j)
>0$. For simplicity, use $Z$-statistic for test and we will reject
each $H_0^j$ if $Z > h$ for some threshold $h$. For level $\alpha
= 0.05$ test, $h=1.64$ (Figure 1). The type I error for the
multiple hypotheses testing would be \bqn\alpha &=& P(\mbox{
reject } H_0|H_0 \mbox{ true })\\
&=& P(\mbox{ reject some } H_0^j|H_0 \mbox{ true })\\
&=& P\Big(\bigcup_{j=1}^m \{Y(x_j) > h\}\Big|\mathbb{E} Y=0\Big).
\eqn Unfortunately, $Y(x_j)$'s are correlated which makes the
computation of type I error almost intractable for random fields
other than Gaussian. Read Keith Worsley's papers at
http://www.math.mcgill.ca/keith. He made random field approach for
multiple comparison very popular in imaging.

 \item {\em i.i.d. simulation.} For simplicity of argument, assume
there is a single image as an observation. For $Z \sim N(0,1)$,
$$P\mbox{-value }= P(Z>1.64) = 0.05.$$ So we expect $5\%$ of
observations to be false positive.
\begin{verbatim}
Y=normrnd(0,1,95,68); imagesc(Y);
colorbar; colormap('bone')
[Yl,Yh]=image_threshold(Y, 1.64);
imagesc(Yh);colorbar; n=95*68
size(find(reshape(Y,n,1)>=1.64),1)/n
>ans =
    0.05
\end{verbatim}


\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.25]{ima17-1.eps}
\includegraphics[scale=0.25]{ima17-2.eps}
\caption{Left: original $N(0,1)$ white noise image. Right: 5$\%$
of pixels are false positives at $\alpha=0.05$ level tests. }
\end{figure}


\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.25]{ima17-3.eps}
\includegraphics[scale=0.25]{ima17-6.eps}
\caption{Left: Gaussian kernel smoothing with 4mm FWHM filter.
FWHM=full width at half maximum. Right: threshold at 1.64 shows no
statistically significant regions. So kernel smoothing can be
considered as an ad-hoc approach for removing false positives.}
\end{figure}


\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.25]{ima17-4.eps}
\includegraphics[scale=0.25]{ima17-5.eps}
\caption{Left: zero mean unit variance Gaussian field constructed
from 12mm fWHM filter. Right: threshold at 1.64 showing
statistically significant clustered regions}
\end{figure}


\item {\em Random field simulation.}

\begin{verbatim}
Y_smoo=gaussblur(Y,4); imagesc(Y_smoo);
colorbar;colormap('bone');
[Yl,Yh]=image_threshold(Y_smoo, 1.64);
imagesc(Yh);colorbar;colormap('bone')

Y_smoo=Y_smoo/std(reshape(Y_smoo,n,1));
>size(find(reshape(Y_smoo,n,1)>=1.64),1)/n
ans =
    0.0424
\end{verbatim}
\end{enumerate}
{\bf Problem 28.} Let $[Y(x_1),\cdots,Y(x_m)]' \sim N(0,\Sigma)$
and $\mathbb{E} [Y(x_i)Y(x_j)] = \sigma_{ij}$. Compute the
probability
$$\alpha(h) = P\big( \mbox{ for some } j,  Y(x_j) > h\big).$$
If you can't solve it analytically, try Monte-Carlo simulation for
many $h$'s and based on the simulation plot $\alpha(h)$.

\end{document}