\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 31\\
Curve Modeling IV.}
\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 {\em Modeling curves.} We show how to simulate
realistic curve. We propose the following model of curve. Let
$X=(x(t),y(t)) \in \mathbb{R}^2$ be the arclength parameterization
of curve $C$.

 $$x(t) = \mu_x(t) + \epsilon_x(t) + w_x(t),$$
 $$y(t) = \mu_y(t) + \epsilon_y(t) + w_y(t).$$
where $\epsilon_x, \epsilon_y$ are independent zero mean Gaussian
random fields and $w_x,w_y$ are independent Gaussian white noises
which is again Gaussian random fields. This is a standard modeling
technique in dependent time series and fMRI signals. So we may let
$\epsilon_x + w_x \to \epsilon_x$ and $\epsilon_y + w_y \to
\epsilon_y$. There is a couple of ways to estimate the mean
coordinate functions $\mu_x$ and $\mu_y$. See principal curve
modeling in Hastie and Tibshirani's The Element of Statistical
Learning (shape analysis version of PCA). For our simple minded
example, we estimate them via iterated-kernel or diffusion
smoothing, i.e.
 $\hat \mu_x = K_{\sigma_x}*x$ and $\hat \mu_y = K_{\sigma_y}*y$
However, if the bandwidth is too large there will be the problem
of shrinkage. See Taubin's seminal paper on the shrinkage effect
of kernel smoothing. Now we take $\hat \mu_x,\hat \mu_y$ to be the
mean coordinate functions. If there are multiple subjects, we need
to perform a curve registration to estimate the mean shape.

\begin{verbatim}
cx=[x(171);x;x(1)]
cy=[y(171);y;y(1)]
for j=1:10
  for i=2:172
    nbr_dist=[dl(i-1) 0 dl(i)];
    W= K(0.2,nbr_dist);
    tx(i-1) = dot(cx((i-1):(i+1)),W);
    ty(i-1) = dot(cy((i-1):(i+1)),W);
  end;
  cx=[tx(171) tx tx(1)];
  cy=[ty(171) ty ty(1)];
end;
\end{verbatim}


\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.24]{ima31-1.eps}
\includegraphics[scale=0.24]{ima31-2.eps}
\includegraphics[scale=0.24]{ima31-3.eps}
\includegraphics[scale=0.24]{ima31-6.eps}
\caption{Top left: iterated-kernel smoothing ($\sigma=0.2$,
$m=10$. Top right: white noise added curve. The white noises have
been exaggerated. Bottom left: Iterated-kernel smoothing of
x-component of noise ($\sigma=0.1$,$m=5$). Bottom right: the
dotted line $(\mu_x,\mu_y)$. the bold line is a simulated curve.}
\end{figure}

\item {\em Simulating curves.} In order to simulate curves, we
need to simulate correlated noise which follows Gaussian random
fields. If we simply add independent Gaussian noise to the
coordinates $w_x(t_i) \sim N(0,\sigma_1^2)$ and $w_y(t_i) \sim
N(0,\sigma_2^2)$ at each point $t_i$ the curve may not be
topologically correct (Figure 1). So we add correlated noise.
$$\epsilon_x(t) = K*w_x(t), \epsilon_y(t) = K*w_y(t).$$
Then simulated curve is given by
$$x(t) = \hat \mu_x + K*w_x, y(t) = \hat \mu_y + K*w_y.$$
Note that $\mathbb{E} x = \hat \mu_x, \mathbb{E} y = \hat \mu_y$.
These guarantee $\hat \mu_x$ and $\hat \mu_y$ are in fact the mean
coordinate functions.

Generalizing this idea, smoothing and simulating random surfaces
follow similarly (Figure 3). Spline-based smoothing and simulation
can be performed as well (Figure 2). But most 2D curve and surface
smoothing methods has the shrinkage problem when the bandwidth is
large.

\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.55]{ima31-4.eps}
\caption{Spline smoothing by Thomas Hoffmann.}
\end{figure}


\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.3]{ima31-5.eps}
\caption{Iterated-kernel smoothing applied to noisy surface.}
\end{figure}

\item {\em Different method.} The above approach uses the fixed
Cartesian coordinate system to add noise components. A different
approach is to add noises along the normal direction. For the
arclength parameterization $X(t)=(x(t),y(t))$, $X'(t)$ and
$X''(t)$ give the unit normal and tangent vectors respectively.
The normalized tangent vector is $N(t) = X''(t)/\|X''(t)\|$. The
Cartesian coordinate representation is
$$N(t) = \frac{(x''(t),y''(t))}{\sqrt{x''(t) + x''(t)}}.$$
If $t$ is interpreted as time, $X,X',X''$ can be viewed as the
position, velocity and acceleration of a particle. Estimating
normal vectors in both noisy curves and surfaces is a nontrivial
problem. This is equivalent to the second order derivative
estimation which is extremely sensitive to noise. In computer
vision literatures, there is a normal vector and curvature
estimation technique called {\em vector voting} or {\em tensor
voting}.

For given mean curve $\mu(t)=(\mu_x(t), \mu_y(t))$, let $N(t) =
(\cos \theta(t), \sin \theta(t))$ be the normal vector of
$\mu(t)$. Add white noise $\|w\| \sim N(0,\sigma_N^2)$ along the
direction of $N(t)$. Then we are adding random vector $w=(w_x,w_y)
=(\|w\|\cos\theta(t),\|w\|\sin \theta(t))$. Note $w_x \sim N(0,
\sigma^2_N\cos^2 \theta(t)), w_y \sim N(0, \sigma_N^2 \cos^2
\theta(t))$. Unlike the previous approach, $w_x$ and $w_y$ are
correlated. Assuming the mean curve gives somewhat smooth normal
vector, we add noise by applying iterated-kernel smoothing on
$\|w\|$. So our stochastic curve model is
$$x(t) = \mu_x(t) + \cos \theta(t)K_{\sigma}*\|w\|(t),$$
$$y(t) = \mu_y(t) + \sin \theta(t)K_{\sigma}*\|w\|(t).$$
Again this also gives $\mathbb{E}x=\mu_x, \mathbb{E}y = \mu_y$.
\end{enumerate}
{\bf Problem 45.} Based on the above method of adding noise in the
normal direction, simulate a coupe of corpus callosum curves. 
\end{document}

\item 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$.
