\pretolerance=1400
\magnification=\magstep1
%\baselineskip=0.77cm

file publ/lzrazky/precip6.tex

nic s tym zatial nerobte len si to precitajte a
poslite mi mi nejaky draft vasej casti

ked to bude v nejakej finalnejsie podobe tak to dam
do wordu alebo html, ale ozaj neviem co so vzorcami
lebo s tym nemam cas sa prplat vo worde

{\bf Precipitation by multivariate RST } 
\bigskip

\noindent
HELENA MITASOVA Department of Geography,
University of Illinois at Urbana-Champaign, Urbana, IL 61801, U.S.A.

\noindent
LUBOS MITAS National Center for Supercomputing Applications,
University of Illinois at Urbana-Champaign, Urbana, IL 61801, U.S.A.

\bigskip
\noindent
{\bf 2.1  Multivariate spline interpolation}

\smallskip
Radial basis functions (multiqudrics, splines) belong
to the most accurate interpolation and approximation
methods for scattered data (REF). 
Because of its accuracy, flexibility, multivariate formulation, 
as well as implementation in GIS, we have selected the 
Regularized Spline with Tension (RST) as one of the methods 
which may be  suitable for interpolation of precipitation 
in a region with significant variability in terrain. 
The RST method has been described by Mitas and Mitasova 1993 
(an earlier version called Completely Regularized Spline)
and Mitasova et al. 1995 (the general d-dimensional formulation with applications
for 2D,3D,and 4D data). Therefore, here we briefly recall only 
the basic principles and focus on new developments and issues 
specific to interpolation of precipitation data.

In general, a spline function $S({\bf x})$ fulfills the
condition of minimizing the deviation from the measured points
and at the same time its smoothness seminorm $I(S)$

$$ 
\sum_{j=1}^N |p^{[j]}-S({\bf x}^{[j]})|^2w_j + w_0I(S)={\sl minimum}\, . 
\eqno(1)
$$

\noindent
where $p^{[j]}$ are the values 
measured at discrete points ${\bf x}^{[j]}=(x_1^{[j]},x_2^{[j]}, \dots ,x_d^{[j]}),
\; j=1,\dots ,N$ within a certain region of a {\sl d}-dimensional space,
$w_j,w_0$ are positive weighting factors and 
$I(S)$ is the measure of smoothness (smooth seminorm) or roughness penalty.
For $w_j/w_0 -> 0$ the  function $S({\bf x})$ passes through the data points
and works as an interpolation function, for  $w_j/w_0 > 0$ the function can
deviate from the data points and  $S({\bf x})$ is an approximation function. 
The general solution of the minimization problem given by Eq. (1)
can be expressed as a sum of two components (Talmi and Gilat 1977)

$$
S({\bf x})=T({\bf x})+\sum_{j=1}^N \lambda_j R({\bf x},{\bf x}^{[j]})
\eqno(2)
$$

\noindent
where $T({\bf x})$ is  a 'trend' function and $R({\bf x},{\bf x}^{[j]})$
is a {\sl  radial basis function} with an explicit form depending
on the choice of the  $I(S)$. 
The smoothness seminorm  $I(S)$ for the RST method has been designed  
to fulfill the following conditions:
\item {a)}  $S({\bf x})$ can be derived in an explicit form;
\item {b)}  $S({\bf x})$ synthesizes in a single function 
desired properties of several previously known splines 
(thin plate spline: Duchon 19XX, spline with tension: Franke 1985, 
Mitas and Mitasova 1988, Hutchinson and xx 1989?, 
regularized spline: Mitas and Mitasova 1988).

\noindent
The seminorm which fulfills these requirements
includes derivatives of {\sl all orders} with weights rapidly 
decreasing with the increasing derivative order (Mitas 1993 unpublished 
manuscript) and  for $d=2$ has the following form:

$$
I^2(S)=\sum_{\bf \alpha} B_{\bf \alpha} \int \int _{\Omega}
\left[ {\partial^{|{\bf \alpha}|}\over
 \partial x_1^{\alpha_1} \partial x_2^{\alpha_2}}
 \, S({\bf x})\right]^2 dx_1 dx_2 
 \eqno(3)
$$

\noindent
where ${\bf \alpha} =(\alpha_1,\alpha_2)$ is a multiindex 
($\alpha_1=0, 1, 2, ...., \, 
\alpha_2=0, 1, 2, ...)$, with $|{\bf \alpha}|= \alpha_1 + \alpha_2$ 
and $\{ B_{\bf \alpha}\} $ are nonnegative constants, in our case 

$$
B_{\bf \alpha}=\cases {\,\quad 0\; ,  \qquad\qquad |{\bf \alpha}| =0 \cr
\quad\cr
 \quad
{|{\bf \alpha}|!\over {\alpha_1! \alpha_2!}}
{1. \over {\varphi ^{2|{\bf \alpha}|}\; (|{\bf \alpha}|-1)!}}\; ,
\quad |{\bf \alpha} |>0
\cr}
\eqno (4)
$$

\noindent
where $\varphi$ is a relative reciprocal weight of particular terms in
the sum (generalized tension)
which provides the control over the influence of derivatives
of certain order on the resulting function.
With this particular choice of coefficients,
 the radial basis functions can be found explicitly
and have the following form for d=2,3,4 (Mitas and Mitasova 1999?
Mitasova et al 1995)

$$
R_2({\bf x},{\bf x}^{[j]}) = R_2(r) =
- \left[ {\rm E_1} (\varrho) + \ln\varrho + C_E \right]
 \, \quad d=2 
\eqno(5)
$$

$$
R_3({\bf x},{\bf x}^{[j]}) = R _ 3(r)=
\sqrt{{\pi \over  \varrho}} \;{\rm erf} \left( \sqrt{\varrho}  \right) -2
\quad d=3
\eqno(6)
$$

$$
R_4({\bf x},{\bf x}^{[j]})=R_4(r)=  {1 -e^{-\varrho} \over \varrho } - 1  \quad d=4 
\eqno(7)
$$

\noindent
where $\varrho = (\varphi r/2)^2$,
 $r^2=\sum_{i=1}^d (x_i-x_i^{[j]})^2$ is the squared distance,
$\varphi$ is the generalized tension, 
$C_E=0.577215...$ is the Euler constant,  ${\rm E}_1(.)$
is the exponential integral function and
erf(.) is the error function (Abramowitz and Stegun 1964).

The coefficients $a_1, \{\lambda _j\}$ are obtained by solving
the following system of linear equations

$$ 
a_1+\sum_{j=1}^N \lambda_j \left[R({\bf x}^{[i]},{\bf x}^{[j]})
+\delta_{ji}w_0/w_j \right] =z^{[i]}\, ,\quad i=1,\ldots,N 
\eqno (8)
$$

$$
\sum_{j=1}^N \lambda _j =0 
\eqno (9)
$$

\noindent
where $w=w_0/w_j???$ is the smoothing parameter.
The computational demands due to the solution of this system
have been (and often still are ) cited as a major
limitation for the use of splines (REF) 
However, the problem has been satisfactorilly resolved
by a segmentation procedure implemented  
in GRASS GIS version of the RST,
which is being routinely used to interpolate data 
sets with hundreds of thousands of points
(Mitasova et al. 1995, Mitas and Mitasova 1999, Hargrove, ...REF). 

\medskip
\noindent
{\bf Anizotropy }
The RST interpolation functions are
(a) invariant to a rotation of  coordinate space
 because the basis functions depend only on distance;
(b) {\sl not} scale invariant, and a change of scale is equivalent
to change in the tension parameter $\varphi$.

The second property can be used for implementation of anisotropy 
by rotating the coordinate system by an angle $\theta$ representing
the direction of anisotropy and then rescaling one axis according
to the anisotropy magnitude.
For $d=2$, the transformation from original
coordinates ${\bf x}=(x_1,x_2)$ to new coordinates ${\bf x'}=(x'_1,x'_2)$ is simply

$$
x'_1 = x_1 cos(\theta) + x_2 sin(\theta)
\eqno (12)
$$
$$
x'_2 = -x_1 sin(\theta) + x_2 cos(\theta) 
\eqno (13)
$$

\noindent
and the distances are then computed as

$$
 r^2 = s (x_1'-x_1'^{[j]})^2+ (x_2'-x_2'^{[j]})^2
\eqno (14)
$$

\noindent
 where $s=s_1/s_2$ is a scaling coefficient expressed here as a ratio
between the scaling in $x'_1$ and $x'_2$ axis respectively. The rotated 
and rescaled distances are then used in the solution of the system of equations (8)
as well as for computation of interpolated values in grid 
points. The direction and scaling of the anisotropy can be found using 
standard geostatistical tools (REF).
The generalization of the coordinate transformation for $d=3$ is
straightforward.

\medskip
{\bf Interpolation with influence of additional variable}
To interpolate precipitation with the influence of topography 
we can use an approach similar to the one proposed by Hutchinson and Bischof
(1983)(OVER) 

Given the $N$ values of precipitation $p^{[j]}$
measured at discrete points ${\bf x}^{[j]}=(x_1^{[j]},x_2^{[j]},x_3^{[j]}),
\; j=1,\dots ,N$ within a certain region of a {\sl 3}-dimensional space
we compute a function $p=F(x_1,x_2)$ (as a grid) 
representing the spatial distribution of precipitation over the terrain
surface $x_3=G(x_1,x_2)$ as

$$
p=F(x_1,x_2)=S(x_1, x_2, cG(x_1,x_2))
\eqno (13)
$$

\noindent
where $c$ is vertical rescaling parameter (vertical tension),
and $S(.)$ is the trivariate RST function.
Equation (13) can be interpreted as an intersection of volume
model of spatial distribution of precipitation
(modeled by trivariate RST) with terrain surface.

For this approach it is not necessary that the correlation between
the elevation and precipitation exists for the entire region - the only
condition is that it is sufficiently represented by data and changes
smoothly (correlation is only local).

\medskip
\noindent
{\bf Control parameters, their importance, tuning and optimization}
One of the advantages of the RST method is its flexibility
(adjustability to various types of phenomena) within the
single radial basis function (as opposed to the rather complex
and often subjective task of selecting the proper semivariogram in kriging).
While most applications require only slight adjustment 
of tension parameters there are applications
(such as precipitation in mountanious regions) when a wider
set of tunable parameters can improve the predictive accuracy.
We investigated the impact of the following parameters on the 
predictive accuracy of RST:
\item {a)} horizontal uniform tension  $\varphi$, 
\item {b)} uniform smoothing $w$,
\item {c)} anisotropy: rotation and scale ($\theta,s$),
\item {d)} vertical tension $c$
\item {e)} level of detail (smoothing) of 3-rd variable

TOTO DALEJ treba napisat podla toho co bude vo vysldekoch -
explain what each of the parameters does

The parameters can be selected empirically,
based on the knowledge of the modelled phenomenon, or automatically
by minimization of the predictive error estimated by
a cross-validation procedure (Mitasova et al 1995). 
(as we show later? minimum CV does not necessarily mean the
best result? daj do discussion/zaverov - consistency and
good results from several measures are important)

Discuss normalization issue and role of tension as a rescaling parameter?.

\medskip
\noindent
{\bf Cross-validation}
Napiste ako ste to robili - toto je text z nasho IJGIS clanku 1995

We have used an ordinary 
cross-validation to find both optimal tension and smoothing.
If we assume that $\{w_j\}$ are the same for all data points then the smoothing
parameter is $w=w_j/w_0$.
Let $S_{N,\varphi,w}^{[k]}$ be the smoothing spline
with tension using all the data points except the k-$th$.
We can take the ability of $S_{N,\varphi,w}^{[k]}$ to
predict the missing data point $z_k$ as a measure
of predictive capability of the spline with the
given tension $\varphi$ and smoothing $w$,
 which can be represented by a mean error $V_0$
 
$$
V_0(\varphi,w)= {1 \over N} \sum _{k=1}^N |S_{N,\varphi,w}^{[k]}
({\bf x}^{[k]})-z^{[k]} |.
 \eqno(15)
$$
\noindent
Then the optimal values of tension $\varphi_O$ and smoothing
$w_O$ are found by  minimizing $V_0(\varphi,w)$. 
The cross-validation error computed for each point can also be used
for the analysis of spatial
distribution of predictive error and location of areas where
this error is high and denser sampling is needed. 

Maybe include some comments on:  many of the geostatistical concepts
 can  be exploited within the spline framework.

Try to fit the RST RBF to variogram? -see whether Roling has done that.


\bigskip
\noindent
{\bf References}
\smallskip

\smallskip\noindent
Abramowitz, M., and  Stegun, I.A., 1964, {\sl Handbook of Mathematical
 Functions} (Dover: New York), 297-300, 228-231.

\smallskip\noindent
Compton, W., and Miller, D., 1994, Extending
environmental modeling into fifth dimension: Multiple concurrent
phenomena. {\sl Proceedings of the II. Int. Conference
on Integration of Environmental Modeling and GIS}, Breckenridge,
Colorado, in press.

\smallskip\noindent
Dietrich, C.R., and Osborne, M.R., 1991, Estimation of covariance 
parameters in kriging via restricted maximum likelihood. {\sl Mathematical
Geology} {\bf 23}, 119-135.

\smallskip\noindent
Franke, R., 1987, Recent advances in the approximation of surfaces
from scattered data. {\sl Topics in Multivariate Approximation},
edited by Chui, Schumaker and Utreras (Academic Press), 79-98.

\smallskip\noindent
{\sl GRASS4.1 Reference Manual}, 1993,
 U.S.Army Corps of Engineers,
Construction Engineering Research Laboratories, Champaign, Illinois, 422-425.

\smallskip\noindent
Hutchinson, M.F. ,1991, The Application of thin plate smoothing splines
to continent-wide data assimilation.
 {\sl Data Assimilation Systems}, edited by Jasper, J.D.,
BMRC Research Report N.27, Bureau of Meteorology, Melbourne, 104-113.

\smallskip\noindent
Hutchinson, M.F., 1994, Interpolating rainfall means - getting the
temporal statistics correct. {\sl Proceedings of the II. Int. Conference
on Integration of Environmental Modeling and GIS}, Breckenridge,
Colorado, in press.

\smallskip\noindent
Hutchinson, M.F., and Bischof, R.J., 1983, A new method for estimating the
spatial distribution of mean seasonal and annual rainfall applied to the
Hunter Valley, New South Wales. {\sl Australian Meteorological magazine}
{\bf 31}, 179-184.

\smallskip\noindent
Hutchinson, M.F. and Gessler, P.E., 1994, Splines -  more than just
a smooth interpolator. {\sl Geoderma} {\bf 62}, 45-67.

\smallskip\noindent
Mitas, L. and Mitasova, H., 1988, General variational approach to
the interpolation problem.
 {\sl Computers and Mathematics with Applications}, {\bf 16}, 983-992.

\smallskip\noindent
Mitas, L. and Mitasova, H., 1994, Multivariate approximation for
scattered data. (to be published).

\smallskip\noindent
Mitasova, H., 1992, 1993, 1994, Surfaces and Modeling.
{\sl GRASSClippings}, {\bf 6}/3, {\bf 7}/1,2,3, p. 16-18, 18-19, 18-19, 25-26.

\smallskip\noindent
Mitasova, H., and  Mitas, L., 1993,
Interpolation by Regularized Spline with Tension: I. Theory  and  implementation.
{\sl Mathematical Geology }, {\bf 25}, 641-655.

\smallskip\noindent Nielson, G.M., Foley, T.A., Hamann, B., and
Lane, D., 1991, Visualizing and Modeling Scattered Multivariate Data.
{\sl IEEE Computer Graphics and Applications} {\bf 11}, 47-55.

\smallskip\noindent
Talmi, A., and Gilat, G., 1977, Method for Smooth
 Approximation of Data. {\sl Journal of Computational Physics} {\bf 23}, 93-123.

\smallskip
\noindent
 {\sl Trends in Nitrogen in the Chesapeake bay, (1984-1990)}, 1992,
 Report CBP/TRS 68/92, Chesapeake Bay Program.

\smallskip\noindent
 Wahba, G., 1990, Spline models for observational data.
{\sl CNMS-NSF Regional Conference series in applied mathematics} {\bf 59}, 
(Philadelphia, Pennsylvania : SIAM).

\bye