\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 35\\
Smoothing periodic functional data}
\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}

\begin{enumerate}
\item Solution to Problem 34. Diffusion smoothing on a unit
circle. Motivation would be when we have periodic functional data
$g(\theta) = g(\theta + L)$. This periodic condition can be
transformed to
$$g\big(\frac{L}{2\pi}\frac{2\pi}{L}\theta  \big)
=g\Big(\frac{L}{2\pi}\big(\frac{2\pi}{L}\theta + 2\pi\big)\Big).$$
So any period $L$ function can be transformed to a periodic
function with period $2\pi$. So without loss of generality, we
consider periodic data of the form $g(\theta) = g(\theta + 2\pi)$
which is viewed as data defined on a unit circle.

Laplacian on $S^1$ is trivially $\Delta = \partial^2_{\theta}$.
This can be seen from the line element $dl = d \theta$ so the
Laplacian on a circle must be the Laplacian in 1D Euclidean space.
We need to solve $\partial_t f = \partial^2_{\theta}f$ with
initial condition $f(\theta,0) = g(\theta)$. The solution to an
isotropic diffusion equation on manifold is given as a heat kernel
convolution by finding the eigenfunctions of Laplacian first. We
solve $\partial_{\theta}^2 H + \lambda H =0$. Trivially two
independent solutions are $H = \cos \sqrt{\lambda} \theta, \sin
\sqrt{\lambda} \theta$ with $\lambda > 0$. To grantee the
periodicity of the eigenfunctions, $\sqrt{\lambda} = j \in
\mathbb{Z}^+$. So we have sequence of eigenfunctions $H_j = \sin
j\theta, \tilde{H}_j=\cos j \theta$. $j=0,1,2,\cdots$. Note that
$\int_0^{2\pi} \sin^2 j \theta \;d\theta= \int_0^{2\pi} \cos^2 j
\theta \;d\theta = \pi$ from symmetry. So we normalize the
eigenfunctions $H_j = \frac{1}{\pi}\sin j\theta, \tilde{H}_j =
\frac{1}{\pi}\cos j\theta$. This gives an orthonormal system of
basis functions on $S^1$.


From Lecture 21, heat kernel expansion, \bqn
K_t(\theta_1,\theta_2) &=& \frac{1}{\pi^2}\sum_{j=0}^{\infty} e^{-
j^2t}[\sin (j\theta_1)\sin (j\theta_2)\\
 &&+ \cos (j\theta_1)\cos
(j\theta_2)]\\
&=& \frac{1}{\pi^2}\sum_{j=0}^{\infty} e^{-j^2t}
\cos[j(\theta_1-\theta_2)].\eqn This makes intuitive sense since
the kernel must be isotropic. Then the solution to the diffusion
smoothing is given by \bqn f(\theta,t) &=& \int_0^{2\pi}
K_t(\theta,\phi)g(\phi)\;d\phi\\
&=& \frac{1}{\pi^2}\sum_{j=0}^{\infty}e^{-j^2t} \int_0^{2\pi}
g(\phi)\cos[j(\theta - \phi)] \;d\phi.\eqn

Obviously it is much easier to discretize the diffusion equation
than implementing the above analytical solution exactly.
\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.25]{ima35-1.eps}
\includegraphics[scale=0.25]{ima35-2.eps}
\caption{Left: White noise $N(0,1)$ and $K_{\sigma}*N(0,1)$.
Right: 1000 zero mean unit variance Gaussian field. We are seeing
the boundary effect of kernel smoothing.}
\end{figure}

\item Solution to Problem 36 (Tulaya Limpiti). Simulate zero mean
isotropic Gaussian field in $[0,1]$ and find the {\em corrected}
P-value numerically. First we generate 5000 Gaussian white noise
$N(0,1)$ apply iterated Gaussian kernel smoothing.
\begin{verbatim}
weight= K(1,[-1 0 1])
>>weight =
    0.2119    0.5761    0.2119
GRF=zeros(5000,99);
for k=1:5000
  w=normrnd(0,1,101,1);
  tem=w;
  for j=1:20
    for i=2:100
      tem(i) = dot(w((i-1):(i+1)),weight);
    end;
    w=tem;
  end;
  w=w(2:100);
  GRF(k,:)=w';
end;
x=[1:99]/99; plot(x,w);
\end{verbatim}
Generating random field in this way cause the boundary effect
effect (Figure 1.). So we need to chop off fields near boundary.
So we will only use {\tt w(10:90)} and rerun the above code.

\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.25]{ima35-3.eps}
\includegraphics[scale=0.25]{ima35-4.eps}
\caption{Left: QQ-plot of  the suprima Gaussian random field.
Right: Monte-Carlo version of $P$-value of the suprima field.}
\end{figure}

\begin{verbatim}
sup=max(GRF,[],2)
pvalue=inline('length(find(sup>=h))/5000')
x=[0:0.01:1.63];
for i=1:length(x)
    P(i)=pvalue(x(i),sup);
end;
plot(x,P,'linewidth',2)
\end{verbatim}

\item {\em Quantile-quantile plot} (QQ-plot) QQ-plot was first
introduced by Wilk and Gnanadesikan (1968). Quantile point $q$ for
random variable $X$ is a point that satisfies $P(X \leq q)
=F_X(q)= p$, where $F_X$ is the cumulative distribution function
(cdf) of $X$. Assuming we can find the inverse cdf,
$q=F_X^{-1}(p)$. There is a way to define the inverse of cdf (see
Casella and Berger's Statistical inference). Then we can define
the quantile-quantile plot of two distributions $X$ and $Y$ to be
a parametric curve $(F_X^{-1}(p), F_Y^{-1}(p)) \in \mathbb{R}^2$
for $p \in [0,1]$.

{\bf Problem 48.} Prove that the QQ plot of $N(\mu_1,\sigma_1^2)$
and $N(\mu_2,\sigma_2^2)$ is a straight line. What is the slope
and the intersection? Bonus problem: Find the QQ plot of two
lognormal distributions $exp(N(\mu_1,\sigma_1^2)),
exp(N(\mu_2,\sigma_2^2))$ analytically.

\end{enumerate}
\end{document}

Consider following random field model on manifolds $Y_i(x) =
\mu_i(x) +\epsilon(x), x \in \partial \Omega$ for subject $j$. We
assume i.i.d. zero mean Gaussian random field $\epsilon_j$ with
isotropic covariance $R(x,y) = \sigma^2f(d(x,y))$ where $d$ is a
geodesic distance and $f$ some function.

Suppose there are other measurement $Z_j \sim N(\mu_Z
,\sigma_Z^2)$. The sample correlation between two measurements
$Y_j$ and $Z_j$ at fixed point $x$ is given by
$$\rho(x) = \frac{\mbox{Cov}(Y(x),Z)}{$$


\item {\em Covariance function estimation.} On manifolds $\partial
\Omega$, we have model
$$Y(x) = \mu(x) + \epsilon(x)$$
where $\epsilon$ is unknown random field. We may estimate $\mu$ by
kernel smoothing so that $\widehat{\mu} = K_{\sigma}*Y$. Now we
estimate the noise with the residual $\widehat{\epsilon} = Y -
K_{\sigma}*Y$.