TITLE: Question on the constraint force of a string
QUESTION [0 upvotes]: I would like to know if the following statement is true or false even if I expect that it is true.
Notation: I will consider a string that has no mass and that can not be extended. Saying that a force has downward trajectory I will mean a force such that goes toward the floor, so this force is not necessarily vertical.
Statement: Consider a vertical string with an extremity fixed on the ceiling. If we attach on the other extremity a point (with or without mass it's the same) and there is a force $\vec F$ acting on this point with downward trajectory, then only the component of $\vec F$ that is parallel to the string "causes" the constraint force $\vec R$ of the string. Furthermore, this $\vec R$ has the same intensity of the parallel component of $\vec F$ that we have just mentioned, but opposite direction.
Moreover, if it is true, can we say that it is an experimental law ?

REPLY [1 votes]: I think you are confusing static and dynamic or are trying to separate the two in a way that isn't helpful and leads to 'over-thinking' of a simple problem. Let's look at the force diagram again:

In the point $P$ a downward (but as you stated not necessarily vertical) force $\vec{F}$ acts on the end of the string. We decompose it into two components: $\vec{F_1}$ which is parallel to the line $OP$ and $\vec{F_2}$ which is perpendicular to the line $OP$.
Assuming the string is strong enough it will provide a reactive force $\vec{R}$, so that $\vec{R}=-\vec{F_1}$.
Now when we look at the forces acting on point $P$ along the axis perpendicular to $OP$ we can see there's a net force acting on $P$, that is the force $\vec{F_2}$. By Newton's second law this means there must now be acceleration along that axis and in the direction of $\vec{F_2}$. It's shown as the red vector $\vec{a}$.
We can't know the magnitude of the vector $\vec{a}$ yet, as no masses have been specified. Some mass must reside in the point $P$ or the string or both (without mass the acceleration would be infinite!)
In addition, as pointed out by philip_0008, as soon as the point $P$ has gained even the tiniest velocity (a consequence of the acceleration it's subject to) a centripetal force $\vec{F_c}$ must act on $P$, along the $OP$ axis and in the direction of $O$. Without it the point $P$ cannot stay on its circular path.
This means that at $t=0$ then $\vec{R}=-\vec{F_1}$, but at $t>0$ then $\vec{F_c}=\vec{R}-\vec{F_1}$.
So, in the absence of forces other than $\vec{F}$ this by definition a dynamic problem: without some counter-acting force motion will always be the result.


TITLE: $p$-adic differentiation
QUESTION [0 upvotes]: $p$-adic norm is not  obtained via inner product. Is it still possible to differentiate $|x|_p$ with respect to $x\in\mathbf{Q}_p$ by thinking of the norm as a function of $x$?

REPLY [0 votes]: I just managed to work this out:

The   Vladimirov derivative $D^{n}f(x)=\frac{1}{Γ_p(−n)}\int\limits_{Q_p}\frac{f(x)−f(y)}{|x−y|^{1+n}_p}d_py$ is indeed defined for $f(x)=|x|_p^{\mu}$, where $\Gamma_p(a)=\frac{1-p^{a-1}}{1-p^{-a}}$ is the $p$-adic gamma function.
To evaluate the derivative  first shift $y$ to $y'=x-y$, then redefine $y''=xy'$. This brings  the $x$-dependent factor out of the integral,
$$D^n|x|^{\mu}_p=\frac{1}{Γ_p(−n)} |x|^{\mu-n}_p\int\limits_{Q_p}\frac{\big(1-|1-y''|^{\mu}_p\big)}{|y''|^{1+n}_p}d_py''.$$
The integral can be evaluated using the usual rules in terms of $p$-adic gamma functions, yielding
$$ D^n |x|^{\mu}_p=\frac{\Gamma_p(1+\mu)}{\Gamma_p(1+\mu-n)}|x|^{\mu-n}_p.$$
Note that this generalises
$$\left(\frac{d}{dx}\right)^nx^{\mu}= \frac{\Gamma(1+\mu)}{\Gamma(1+\mu-n)}x^{\mu-n}$$
on the real line.
Thank you for your responses.


text:
\begin{document}

\begin{center}
{\Large Valid model-free spatial prediction}
\\\vspace{6pt}
{\large Huiying Mao\footnote[1]{The Statistical and Applied Mathematical Sciences Institute}$^,$\footnote[2]{North Carolina State University} \; Ryan Martin$^2$ \; and Brian J Reich$^2$}\\
\today
\end{center}

\begin{abstract}\begin{singlespace}
\noindent Predicting the response at an unobserved location is a fundamental problem in spatial statistics.  Given the difficulty in modeling spatial dependence, especially in non-stationary cases, model-based prediction intervals are at risk of misspecification bias that can negatively affect their validity. Here we present a new approach for model-free spatial prediction based on the {\em conformal prediction} machinery.  Our key observation is that spatial data can be treated as exactly or approximately exchangeable in a wide range of settings.  For example, when the spatial locations are deterministic, we prove that the response values are, in a certain sense, locally approximately exchangeable for a broad class of spatial processes, and we develop a local spatial conformal prediction algorithm that yields valid prediction intervals without model assumptions. Numerical examples with both real and simulated data confirm that the proposed conformal prediction intervals are valid and generally more efficient than existing model-based procedures across a range of non-stationary and non-Gaussian settings.
\vspace{12pt}\\
{\bf Key words:} Conformal prediction; Gaussian process; Kriging; non-stationary; plausibility.  \end{singlespace}\end{abstract}

\section{Introduction}\label{s:intro}

Providing valid predictions of the response at an unobserved location is a fundamental problem in spatial statistics.  For example, epidemiologists may wish to extrapolate air pollution concentrations from a network of stationary monitors to the residential locations of the study participants.  There are a number of challenges one faces in carrying out valid prediction at a new spatial location, but one of the most pressing is that existing methods are model-based, so the reliability of the predictions depends crucially on the soundness of the posited model. For example, prediction intervals based on Kriging---see \citet{cressie1992statistics} and Section~\ref{s:back:spatial}---often rely on normality, and stationarity is often assumed to facilitate estimating the spatial covariance function required for Kriging. It is now common to perform geostatistical analysis for massive datasets collected over a vast and diverse spatial domain \citep[e.g.,][]{heaton2019case}.  For complex processes observed over a large domain, the normality and stationarity assumptions can be questionable. Failing to account for nonstationarity can affect prediction accuracy, but typically has a larger effect on uncertainty quantification such as prediction intervals \citep{fuglstad2015does}.  While there are now many methods available for dealing with nonstationary \citep[see][for a recent review]{risser2016nonstationary} and non-Gaussianity  \citep[e.g.,][]{gelfand2005bayesian,duan2007generalized,reich2007multivariate,rodriguez2011nonparametric}, these typically involve heavy computations.  This exacerbates the already imposing computational challenges posed by massive datasets. Further, fitting the entire stochastic process may be unnecessary if only prediction intervals are desired.  Therefore, having a method that provides provably valid prediction intervals without requiring specification of a statistical model---and, hence, not inheriting the risk of model misspecification biases---would be very useful to practitioners.  Developing such a method is the goal of this paper.

In recent years, the use of machine learning techniques in statistics has become increasingly more common.  While there are numerous examples of this phenomenon, the one most relevant here is {\em conformal prediction}. This method originated in \citet{vovk2005algorithmic} and the references therein \citep[see, also,][]{shafer2008tutorial}, but has appeared frequently in the recent statistics literature \citep[e.g.,][]{lei2014distribution, lei2018distribution, guan2019conformal, romano2019conformalized, conformal.shift}. What makes this method especially attractive is that it provides provably valid prediction intervals without specification of a statistical model.  More precisely, the conformal prediction intervals achieve the nominal frequentist prediction coverage probability, uniformly over all data distributions; see  Section~\ref{s:back_conformal}.  The crucial assumption behind the validity of conformal prediction is that the data are exchangeable.

Whether it is reasonable to assume exchangeability in a spatial application depends on how the data are sampled.  On the one hand, if the locations are randomly sampled in the spatial domain, then exchangeability holds automatically; see Lemma~\ref{lem:iid.locations}.  Therefore, the standard conformal prediction framework can be used basically off the shelf.  On the other hand, if the locations are fixed in the spatial domain, then exchangeability does not hold in general.  But we show that, for a wide range of spatial processes, the response variables at tightly concentrated locations are approximately exchangeable; see Theorem~\ref{thm:local.exchangeable}.  Therefore, a version of the basic conformal prediction method applied to these tightly concentrated observations ought to be approximately valid.

Using this insight about the connection between exchangeability and the sampling design, we propose two related spatial conformal prediction methods.  The first, a so-called {\em global spatial conformal prediction} (GSCP) method, described in Section~\ref{s:GSCP}, is designed specifically for cases where the spatial locations are sampled at random.  In particular, this global method produces a prediction interval which is marginally valid, i.e., valid on average with respect to the distribution of the target location at which prediction is desired.  The second, a {\em local spatial conformal prediction} (LSCP) method, described in Section~\ref{s:LSCP}, is designed specifically for the case when the spatial locations are fixed.  Since our goal is to proceed without strong assumptions about the spatial dependence structure, it is only possible to establish approximate or local exchangeability.  Therefore, the proposed local spatial conformal prediction method can only provide approximately valid predictions; see Theorem~\ref{thm:valid.local}.  But our goal in the fixed-location case resembles the ``conditional validity'' target in the conformal prediction literature \citep[e.g.,][]{barber2019limits} so, given the impossibility theorems in the latter context, approximate validity is all that can be expected.

For both the global and local formulations, our proposed method is model-free, robust, and computationally feasible for large datasets. In Sections~\ref{s:sim} and \ref{s:app}, we show using real and simulated data that the proposed methods outperform both standard global Kriging and local approximate Gaussian process regression \citep[laGP;][]{gramacy2015local} for non-stationary and non-Gaussian data.  In addition to be useful for spatial applications, it is also a major advance in conformal prediction to the case of dependent data, and establishes the conditions on the spatial sampling design and data-generating mechanism that ensure validity---or at least approximate validity---of the conformal prediction intervals.

\section{Background} \label{s:back}

\subsection{Spatial prediction methods} \label{s:back:spatial}

Let $Y_i = Y(s_i)\in \RR$ be the observation at spatial location $s_i\in\calD\subseteq \RR^2$.  The data consist of $n$ observations $Z_i=(s_i,Y_i)$.  Geostatistical analysis often assumes that the data follow a Gaussian process with stationary and isotropic covariance function.  That is, $Y_i = \theta_i + \varepsilon_i$, where $\theta_i$ is a mean-zero Gaussian process with covariance $\cov(\theta_i,\theta_j) = \sigma^2 \rho(d_{ij})$, a function of the distance $d_{ij}$ between locations $s_i$ and $s_j$, and $\varepsilon_i\iid\mbox{Normal}(0,\tau^2)$ is error.  A common example is the $\Matern$ correlation function $$\rho(d; \phi, \kappa)=2^{1-\kappa} \Gamma(\kappa)^{-1} (d \phi^{-1} \sqrt{2\kappa})^\kappa K_\kappa(d \phi^{-1} \sqrt{2\kappa}), \quad d \geq 0,$$ where $\phi$ is the correlation range, $\kappa$ is the smoothness parameter, $K_\kappa$ is the modified Bessel function of the second kind.  Covariates can easily be added to the mean of $Y_i$ but for simplicity we assume the process has mean zero.  The main assumptions of this model are that the data are Gaussian and the covariance function is stationary and isotropic, i.e., it is a function only of the distance between spatial locations and is thus the same across the spatial domain.

To make prediction inference at the new location $s_{n+1}$, we study the random variable $Y_{n+1} = Y(s_{n+1})$ together with the other $Y_i$'s.  The full data set, including response to be predicted, $Y_{n+1}$, is denoted $Z^{n+1} = \{Z_1,...,Z_n, Z_{n+1}\}$ and the collection of data excluding $Y_i$ (but including $s_i$) is denoted $Z_{(i)}^{n+1}$. Note that $Z_{(n+1)}^{n+1}$ represent the full observed data set $\{Z_1, \ldots, Z_n\}$ and $s_{n+1}$.

The joint distribution of $Y^{n+1}=(Y_1,...,Y_n, Y_{n+1})$, given $(s_1,...,s_n,s_{n+1})$ and the spatial covariance parameters $\Theta=\{\sigma^2,\tau^2,\phi, \kappa\}$, is $\mbox{Normal}\big(0^{n+1},\Sigma(\Theta)\big)$, where $0^{n+1}$ is a vector of zeros with length $n+1$ and $\Sigma(\Theta)$ is a $(n+1)\times(n+1)$ covariance matrix with $(i,j)$ element being $\sigma^2\rho(d_{ij}; \phi, \kappa) + I(i=j)\tau^2$.
Let ${\hat \Theta}$ be an estimate of $\Theta$ and ${\hat Q} = \Sigma({\hat \Theta})^{-1} = \{{\hat q}_{ij}\}$.  The Kriging prediction and variance of $Y_i$, given $\Theta  = {\hat \Theta}$ and $Z_{(i)}^{n+1}$, are
\begin{equation}
\label{eq:kriging}
{\hat Y}_i = \E(Y_i \mid Z_{(i)}^{n+1})=  -\hat q_{ii}^{-1} \sum_{j\ne i}{\hat q}_{ij}Y_j \quad \text{and} \quad
\hat v_i= \var(Y_i \mid Z_{(i)}^{n+1}) = {\hat q}_{ii}^{-1}, \quad i=1,\ldots,n+1.
\end{equation}
Assuming normality and stationarity gives the predictive distribution $(Y_i \mid Z_{(i)}^{n+1}) \sim\mbox{Normal}({\hat Y}_i,\hat v_i)$ and standardized residuals
\begin{equation}
\label{eq:residual}
e_i = \hat v_i^{-1/2} (Y_i-{\hat Y}_i), \quad i=1,\ldots,n+1.
\end{equation}
In particular, a $100(1-\alpha)\%$ prediction interval for $Y_{n+1}$ is $\hat Y_{n+1} \pm z_{\alpha/2} \hat v_{n+1}^{1/2}$, where $z_{\alpha/2}$ is the upper $\alpha/2$ quantile of a standard normal distribution.

\subsection{Conformal prediction} \label{s:back_conformal}

In this section we take a step back from spatial data analysis and review conformal prediction for non-spatial problems.  Suppose we have a data sequence $Y_1, \ldots, Y_n, Y_{n+1}, \ldots$ assumed to be exchangeable with (joint) distribution $\prob$.  No assumptions about $\prob$ are made here, beyond that it is exchangeable. We observe values $y^n = (y_1,\ldots,y_n)$ of $Y^n = (Y_1,\ldots,Y_n)$ and the goal is to predict $Y_{n+1}$.  More specifically, we seek a procedure that returns, for any $\alpha \in (0,1)$, a prediction interval  $\Gamma^\alpha(y^n)\subseteq \RR$ that is {\em valid} in the sense that
\begin{equation}
\label{eq:valid}
\inf \prob^{n+1}\{Y_{n+1}\in\Gamma^\alpha(Y^n)\}
\geq 1-\alpha, \quad \text{for all $n$},
\end{equation}
where $\prob^{n+1}$ represents the distribution of $(Y_1,\ldots,Y_n,Y_{n+1})$ under $\prob$, and the infimum is over all exchangeable distributions $\prob$.  That we require attaining the nominal frequentist coverage for all $\prob$ rules out the use of model-based procedures, such as likelihood or Bayesian methods.

The original conformal prediction method proceeds as follows.  Define a {\em non-conformity measure} $\Delta(B,z)$, a function that takes two arguments: the first is a ``bag'' $B$ that consists of a finite collection of data points; the second is a single data point $z$.  Then $\Delta(B,z)$ measures how closely $z$ represents the data points in bag $B$.  For example, if $\pi$ is a prediction rule and $d$ is some measure of distance, then we might take $\Delta(B,z) = d(\pi_B, z)$, the distance between $z$ and the value $\pi_B$ returned by the prediction rule $\pi$ applied to $B$.  The choice of non-conformity measure $\Delta$ depends on the context and is at the discretion of the data analyst; often there is a natural choice, as in our spatial prediction setting described in Section~\ref{s:GSCP}.  The only technical requirement is symmetry in the first argument; that is, shuffling the data in bag $B$ does not change the value of $\Delta(B,z)$.

Given the non-conformity measure $\Delta$, the next step is to transform the data in a suitable way via $\Delta$.  Specifically, augment the observed data $y^n = (y_1, \cdots, y_n)$ with a provisional value $y_{n+1}$ of $Y_{n+1}$; this $y_{n+1}$ value is generic and free to vary.  Define
\[ \delta_i = \Delta(\{y_1,\ldots,y_{n+1}\} \setminus \{y_i\}, y_i), \quad i=1,\ldots,n,n+1. \]
Note that $\delta_{n+1}$ is special because it compares the actual observed data with this provisional value of the unobserved future observation.  Next, compute the {\em plausibility} of $y_{n+1}$ as a value for $Y_{n+1}$ according to the formula
\begin{equation}\label{eq:plausibility}
p(y_{n+1} \mid y^n) = \frac{1}{n+1} \sum_{i=1}^{n+1} 1\{\delta_i \geq \delta_{n+1}\},
\end{equation}
where $1\{A\}$ denotes the indicator of event $A$.  Note that this process can be carried out for any provisional value $y_{n+1}$, so the result is actually a mapping $\tilde y \mapsto p(\tilde y \mid y^n)$ which we will refer to as the {\em plausibility contour} returned by the conformal algorithm.  This function can be plotted (see Figure \ref{f:plausibility}) to visualize the uncertainty about $Y_{n+1}$ based on the given data $y^n$, the choice of non-conformity measure, etc.  Moreover, a prediction set $\Gamma^\alpha(y^n)$ can be obtained as the upper level set of the plausibility contour,
\[ \Gamma^\alpha(y^n) = \{\tilde y: p(\tilde y \mid y^n) \geq t_n(\alpha)\}, \]
with threshold $t_n(\alpha) = (n+1)^{-1}\lfloor(n+1)\alpha\rfloor$, which approximately equals $\alpha$ when $n$ is large.

For a quick explanation of why this works, i.e., why \eqref{eq:valid} holds, first note that the event $\{Y_{n+1} \in \Gamma^\alpha(Y^n)\}$ corresponds to $\{p(Y_{n+1} \mid Y^n) \geq \alpha\}$.  So to calculate the probability in \eqref{eq:valid}, we need to know something about the distribution of the random variable $p(Y_{n+1} \mid Y^n)$.  Towards this, the key point is that the provisional value being used in this latter expression is exactly the $Y_{n+1}$ being predicted.  Then exchangeability of the original data sequence and the permutation invariance of $\Delta$ implies that $\delta_1,\ldots,\delta_n,\delta_{n+1}$ are exchangeable too.  And since $p(Y_{n+1} \mid Y^n)$ is proportional to the rank of $\delta_{n+1}$, the claim in \eqref{eq:valid} follows from the elementary result that the rank of any individual in a set of exchangeable random variables has a discrete uniform distribution.

\begin{figure}[t]
\centering

\includegraphics[width=0.7\textwidth]{plausibility_demo.pdf}
\caption{Plot of a plausibility contour and a visualization of two level sets corresponding to $\alpha=0.05$ (red dash line) and $\alpha=0.2$ (blue dotted line).}
\label{f:plausibility}
\end{figure}

Conformal prediction is somewhat mysterious---it is not obviously part of classical or Bayesian statistics.  Some authors  \citep[e.g.,][]{lei2018distribution} refer to $p(y_{n+1} \mid y^n)$ as a p-value corresponding to a ``test'' of the hypothesis that $Y_{n+1}=y_{n+1}$.  This interpretation is potentially confusing because we are accustomed with testing hypotheses about fixed unknowns, not about random variables.  The plausibility notion we are adopting here dates back to \citet{shafer1976mathematical}, but has reemerged recently in a new framework for inference and prediction, namely, {\em inferential models} \citep{martin2013inferential, imbook, martin2016pred, martin2019false}.  More recently, \citet{cella2020valid} have shown that the mapping in \eqref{eq:plausibility} is precisely the plausibility contour resulting from a suitable inferential model construction.  Therefore, the validity property \eqref{eq:valid} achieved by conformal prediction is a special case of the general inferential model validity theory.  For our purposes, we will not need any details from the inferential model framework, but we believe that this characterization is helpful for understanding what conformal prediction is and why it works.

\section{Global spatial conformal prediction} \label{s:GSCP}

\subsection{Algorithm GSCP}

The first approach we consider is a direct application of the original conformal algorithm to spatial prediction but with spatial dependence encoded in the non-conformity measure.  In contrast to the local algorithm presented in Section~\ref{s:LSCP}, this method equally weights the non-conformity across all spatial locations in the plausibility contour evaluation.  Therefore, we refer to this as global spatial conformal prediction, or GSCP for short.

From Section~\ref{s:back:spatial}, recall the notation $Z_i=(s_i, Y_i)$ and $Z^{n+1} = \{Z_1,\ldots,Z_{n+1}\}$; also, $Z_{(i)}^{n+1}$ denotes $Z^{n+1}$ with $Y_i$ excluded but $s_i$ included.  Now define the non-conformity measure for the GSCP algorithm as
\begin{equation}\label{e:A_stat}
 \delta_i =  \Delta(Z_{(i)}^{n+1},Y_i) = d(\hat Y_i, Y_i), \quad i=1,\ldots,n+1,
\end{equation}
where $\hat Y_i$ is the point prediction at spatial location $s_i$, and $d$ measures the distance between $\hat Y_i$ and $Y_i$.  Then we define a plausibility contour exactly like in Section~\ref{s:back_conformal}, with obvious notational changes.  That is, for a provisional value $y_{n+1}$ of $Y_{n+1}$, we have
\begin{equation}\label{eq:GSCP:plausibility}
 p(y_{n+1} \mid Z_{(n+1)}^{n+1}) = \frac{1}{n+1} \sum_{i=1}^n 1\{\delta_i \geq \delta_{n+1}\}.
\end{equation}
And the $100(1-\alpha)$\% conformal prediction interval for $Y_{n+1}$, denoted by $\Gamma^\alpha(Z_{(n+1)}^{n+1})$, is just an upper level set of the plausibility contour, consisting of all those provisional $y_{n+1}$ with plausibility exceeding the $t_n(\alpha)$ threshold.  Algorithm \ref{ag:gscp} states the specific steps of Algorithm GSCP.

Any reasonable choice of point prediction procedure can serve the purpose here, including inverse distance weighting (IDW) predictions \citep{henley2012nonparametric}, Kriging predictions, etc. There is also flexibility in choosing the  distance measure.  The results are invariant to monotone transformations of the distance measure so that, e.g., $d(\hat Y_i, Y_i) = (\hat Y_i - Y_i)^2$ and $d(\hat Y_i, Y_i) = |\hat Y_i - Y_i|$ will give the same plausibility contour, and therefore role of the distance measure is merely to rank prediction errors.  In our numerical results presented below, we use the standardized Kriging residuals $\delta_i = e_i^2$ from \eqref{eq:residual} as the non-conformity measure, but the results are similar for other measures including unstandardized squared residuals.  Appendix~\ref{s:A1} shows how to quickly compute the plausibility contour and prediction interval by exploiting the inherent quadratic structure of the Kriging-based non-conformity measure.

\begin{algorithm}[t]
\SetAlgoLined
\KwIn{observations $z_i = (s_i, y_i)$, $i=1,\ldots,n$; prediction location $s_{n+1}$; non-conformity measure $\Delta$; significance level $\alpha$; and a fine grid of candidate $y_{n+1}$ values}
\KwOut{$100(1-\alpha)\%$ prediction interval, $\Gamma^\alpha = \Gamma^\alpha(z_{(n+1)}^{n+1})$, for $Y_{n+1} = Y(s_{n+1})$}

\For{$y_{n+1}$ in the specified grid}{

\For{$i=1,\ldots,n+1$}{
define $z_{(i)}^{n+1}$ by removing $y_i$ from $z^{n+1}$\;

$\delta_i \gets \Delta(z_{(i)}^{n+1}, y_i)$\;

}

compute plausibility for $y_{n+1}$ as $p(y_{n+1} \mid z_{(n+1)}^{n+1})$ in \eqref{eq:GSCP:plausibility}\;
include $y_{n+1}$ in $\Gamma^\alpha$ if $p(y_{n+1} \mid z_{(n+1)}^{n+1}) \geq t_n(\alpha)$\;
}

\Return $\Gamma^\alpha$.

\caption{Global spatial conformal prediction (GSCP).} \label{ag:gscp}
\end{algorithm}

\subsection{Theoretical validity of GSCP}
\label{s:theory_global}

To start, consider randomly sampled spatial locations, i.e., $S_1,\ldots,S_n$ independent and identically distributed (iid) in $\calD$.  Then write $Z_i = (S_i, Y(S_i))$ for $i=1,2,\ldots$ as before.  Note that the distribution, say $\prob$, of the $Z_1,Z_2,\ldots$ process is determined by that of the $S_i$'s and of the spatial process $Y$.  Given that exchangeability is crucial to the validity of the conformal prediction method, and the non-exchangeability of the distribution of $Y$ given the $S_i$'s, it is natural to ask if the $Z_i$ process defined above is exchangeable.  The following lemma answers this question in the affirmative.

\begin{lem}
\label{lem:iid.locations}
If the spatial locations $S_1,S_2,\ldots$ are iid, then $Z_1,Z_2,\ldots$, with $Z_i = (S_i, Y(S_i))$, is an exchangeable sequence.
\end{lem}

\begin{proof}
Conditioning on $Y$, it is easy to see that $Z_i = (S_i, Y(S_i))$ are iid.  Then the claim follows since conditionally iid is equivalent to exchangeability \citep[e.g.,][]{hewitt1955symmetric}.
\end{proof}

Having established that randomly sampled spatial locations implies that the data are exchangeable, a validity property for GSCP follows immediately from the general theory in, e.g., \citet{shafer2008tutorial}.  The following theorem gives a precise statement.  To highlight the prediction interval's dependence on the location $S_{n+1}$, we use the notation $\Gamma^\alpha(Z^n; S_{n+1})$.

\begin{thm}
\label{thm:valid.global}
Let $Y$ be a stochastic process over $\calD$ and suppose that $S_1,S_2,\ldots$ in $\calD$ are iid.  Define $Z_i = (S_i, Y(S_i))$ for $i=1,2,\ldots$.  Then the prediction interval for $Y(S_{n+1})$ returned by Algorithm GSCP is valid in the sense that
\begin{equation}
\label{e:mv}
\inf \prob^{n+1}\{\Gamma^\alpha(Z^n; S_{n+1}) \ni Y(S_{n+1})\} \geq 1-\alpha, \quad \text{for all $\alpha \in (0,1)$ and $n$},
\end{equation}
where $\prob^{n+1}$ is the joint distribution of $Z_1,\ldots,Z_n,Z_{n+1}$ under $\prob$, and the infimum is over all such exchangeable laws.
\end{thm}

Note that Theorem~\ref{thm:valid.global} makes no assumptions about the spatial process $Y$, so the conclusions hold for non-Gaussian and non-stationary processes.  However, Theorem~\ref{thm:valid.global} only gives a validity result which is marginal in the sense that it accurately predicts the response $Y(S_{n+1})$ at a randomly sampled spatial location $S_{n+1}$.  It would also be interesting to ask if, still under the randomly sampled spatial location regime, the conformal prediction is conditionally valid, i.e., given $S_{n+1}=s^\star$, where $s^\star$ is a fixed spatial location.  There are negative results in the literature \citep[e.g.,][]{lei2014distribution} which state that strong conditional validity---for all $\prob$ and almost all targets $s^\star$---is impossible with conformal prediction.  A considerable amount of effort has been expended recently trying to achieve ``approximate'' conditional validity in some sense; see, e.g., \citet{lei2014distribution}, \citet{conformal.shift}, and \citet{barber2019limits}.

Surprisingly, there is at least one scenario in which a strong conditional validity result can be achieved in our context.  It turns out that one of the obstacles that prevents a conditional validity result in the existing literature is an ``edge effect.''  That is, conditional validity is typically achieved at targets in the middle of the domain, but fails at targets in the extremes; see Figure~1(b) in \citet{lei2014distribution}.  So if it were possible to eliminate the edge effect---even if in a trivial way, by eliminating the edge itself---then there is hope for establishing a conditional validity result.  In our spatial context, but perhaps not in other cases, it may not be unreasonable to assume that the spatial locations are sampled iid from a uniform distribution on a sphere.  Since the sphere has no edges and a uniform distribution has no extremes, there is no ``edge effect,'' and, hence, conditional validity can be established.  Some additional structure in the $Y$ process is also needed here, in particular, it should be isotropic in the sense that the correlation structure only depends on the distance between spatial locations.

\begin{prop}
\label{prop:sphere}
Let $Y$ be an isotropic stationary process over the sphere $\calD = \{s \in \RR^3: \|s\|=1\}$, and suppose that the spatial locations $S_1,\ldots,S_n, S_{n+1}$ are independent and uniformly distributed on $\calD$.  Define $Z_i = (S_i, Y(S_i))$ for $i=1,\ldots,n,n+1$.  Then the prediction interval for $Y(S_{n+1})$ returned by the GSCP algorithm is conditionally valid in the sense that
\[ \inf \prob^{n+1}\{\Gamma^\alpha(Z^n; S_{n+1}) \ni Y(S_{n+1}) \mid S_{n+1} = s^\star\} \geq 1-\alpha \quad \text{for all $s^\star \in \calD$, $\alpha \in (0,1)$, and $n$}, \]
where $\prob^{n+1}$ is the joint distribution of $Z_1,\ldots,Z_n,Z_{n+1}$ under $\prob$, and the infimum is over all such exchangeable laws with uniform marginals for the spatial locations.
\end{prop}

\begin{proof}
See Appendix~\ref{proof:sphere}.
\end{proof}

To our knowledge, the above proposition gives the first conditional validity result for conformal prediction in the literature.  We were able to side-step the usual difficulties by taking advantage of the special structure of the sphere and the uniform distribution thereon. However, these structures, along with the assumed stationarity of $Y$, limit our flexibility.  Fortunately, we can regain that flexibility with the local algorithm to be described in Section~\ref{s:LSCP}.

\section{Local spatial conformal prediction} \label{s:LSCP}

\subsection{Local exchangeability}

Let $\calD$ be the spatial domain, which we take to be $[0,1]^2$ for our analysis here.  Suppose that the $Y$ process defined over $\calD$ can be decomposed as
\begin{equation}
\label{eq:decomp}
Y(s) = \psi\bigl( X(s), E(s) \bigr), \quad s \in \calD,
\end{equation}
where $X$ and $E$ are independent, $X$ is a continuous spatial process, $E$ is a non-spatial process, and $\psi$ is a deterministic, continuous, and real-valued function, such as $\psi(x,e) = x + e$ or $\psi(x,e) = xe$.

Define the sequence of {\em localized processes} around $s^\star$ as
\begin{equation}
\label{eq:local}
\tilde Y_r(u) = Y(s^\star + r u), \quad u \in \calU = \{u \in \RR^2: \|u\| \leq 1\},
\end{equation}
indexed by the unit disk $\calU$ and the radius $r > 0$.  Similarly, write $\tilde X_r$ and $\tilde E_r$ for the localized versions of the $X$ and $E$ processes.  We will require the following:
\begin{itemize}
\item $X$ is {\em $L_2$-continuous} at $s^\star$, i.e., $\E\{X(s) - X(s^\star)\}^2 \to 0$ as $s \to s^\star$; and
\item $E$ is {\em locally iid} at $s^\star$, i.e., $\tilde E_r$ converges in distribution to an iid process as $r \to 0$.
\end{itemize}
This general formulation covers a wide range of common models for a continuous response, including the additive model in Section~\ref{s:back:spatial}, certain generalized spatial linear models \cite{diggle1998model}, spatial copula models \citep{krupskii2019copula}, and max-stable processes \citep{reich2012hierarchical}.  For example, in generalized spatial linear model, with suitable spatial process $X$ and Gaussian white noise $E$, take
\begin{equation} \label{eq:glm}
  \psi(x,e) = H_{g(x)}^{-1}\bigl( \Phi^{-1}(e) \bigr),
\end{equation}
where $H_\xi$ is the distribution function for the exponential family with natural parameter $\xi$, $g$ is the chosen link function, and $\Phi$ is the standard normal distribution function.  Of course, to achieve the continuity requirement for $\psi$, which is essential, the exponential family must have a density with respect to Lebesgue measure.

These assumptions yield a certain kind of {\em local exchangeability} which can then be used to show that our local spatial conformal prediction achieves a desired validity property.  We first establish this local exchangeability result, which may be of independent interest.

\begin{thm}
\label{thm:local.exchangeable}
Suppose that $Y$ can be decomposed as in \eqref{eq:decomp}, where $X$ is $L_2$-continuous at $s^\star$ , $E$ is locally iid at $s^\star$, and $X$ and $E$ are independent. Then the localized process $\tilde Y_r$ in \eqref{eq:local} converges in distribution as $r \to 0$, and the limit is an exchangeable process in the sense that its finite-dimensional distributions are exchangeable.
\end{thm}

\begin{proof}
See Appendix~\ref{proof:local.exchangeable}.
\end{proof}

Theorem~\ref{thm:local.exchangeable} says that, if the spatial locations are situated sufficiently close to $s^\star$, then the joint distribution of the corresponding response values is approximately exchangeable and, therefore, a local spatial conformal prediction approach---LSCP for short---that is based on only those data points sufficiently close to $s^\star$ ought to be valid.  The next section formally states the LSCP algorithm and establishes its theoretical validity properties.

\subsection{LSCP algorithm and its theoretical validity}

Recall that $\calD = [0,1]^2$ is the unit square.  Take a uniform grid $\{N^{-1}, 2 N^{-1}, \ldots, 1\}^2 \subset \calD$ of spatial locations and denote these as $s_i$ for $i=1,\ldots,n$, with $n=N^2$.  We assume that the $Y$ process can be observed on this deterministic grid of spatial locations.  As is common in the spatial statistics literature \citep[e.g.,][]{stein1990comparison,cressie1992statistics}, we adopt an {\em infill asymptotic regime} in which $N \to \infty$.  That is, the region $\calD$ remains fixed while the number of observations in that region increases, filling the space.  The relevant point for our analysis is that under this regime the number of locations at which $Y$ is observed in any neighborhood of $s^\star$ will go to infinity.  Such a regime is natural in cases like ours where without structural assumptions about $Y$, e.g., stationarity, it is simply not possible to learn the local features of a process at $s^\star$ if data are not concentrated in a neighborhood around $s^\star$.

Take a neighborhood around $s^\star$ of radius $r$.  The number of spatial grid points in that neighborhood is roughly $\pi n r^2$, so we need to choose $r$ so that this number is sufficiently large.  In particular, to ensure that the neighborhood contains at least $m$ grid points, where $m > 0$ is some fixed integer, we roughly need
\begin{equation}
\label{eq:radius}
r \propto \Bigl( \frac{m}{n} \Bigr)^{1/2} = \frac{m^{1/2}}{N} \to 0 \quad \text{as $N \to \infty$}.
\end{equation}
It follows from Theorem~\ref{thm:local.exchangeable} that the joint distribution of the response $Y$ at these $m$-many spatial locations around $s^\star$ (corresponding to $m$-many vectors in $\calU$) would be approximately exchangeable and, consequently, a conformal prediction algorithm that creates non-conformity scores using only these $m$ observations would be valid for predicting $Y(s^\star)$.

More precisely, form $Z_i = (s_i, Y_i)$, for $i=1,\ldots,m$, based on only those $m$ locations closest to $s^\star$.  Set $s_{m+1}=s^\star$ and take $y_{m+1}$ as a provisional value for $Y_{m+1}=Y(s^\star)$.  Using a suitable prediction procedure, like that based on Kriging discussed above, define
\[ \delta_i = \Delta(Z_{(i)}^{m+1}, Y_i), \quad i=1,\ldots,m+1, \]
and the corresponding plausibility contour
\begin{equation} \label{eq:LSCP:plausibility}
    p(y_{m+1} \mid Z_{(m+1)}^{m+1}) = \frac{1}{m+1} \sum_{i=1}^{m+1} 1\{\delta_i \geq \delta_{m+1}\}.
\end{equation}
Specific details are presented in Algorithm~\ref{ag:lscp}.

\begin{algorithm}[t]
\SetAlgoLined
\KwIn{observations $z_i = (s_i, y_i), i=1,\ldots,n$; predict location $s^\star$; non-conformity measure $\Delta$; significance level $\alpha$; and a fine grid of candidate $y^\star$ values}
\parameter{number of neighbors to consider $m \leq n$}
\KwOut{$(1-\alpha)100\%$ prediction interval, $\Gamma^\alpha$, for $Y(s^\star)$: }

form $z_i, i = 1, \ldots, m$, based on $m$ locations closest to $s^\star$\;

$s_{m+1} \gets s^\star$\;

\For{$y_{m+1}$ in the specified grid}{

\For{i = 1 to m+1}{
define $z^{m+1}_{(i)}$ by removing $y_i$ from $z^{m+1}$\;
$\delta_i \gets d(z^{m+1}_{(i)}, y_i)$\;}
compute plausibility for $y_{m+1}$ as  $p(y_{m+1} \mid z_{(m+1)}^{m+1})$ in \eqref{eq:LSCP:plausibility}\;

include $y_{m+1}$ in $\Gamma^\alpha$ if $p(y_{m+1} \mid z_{(m+1)}^{m+1}) \geq t_m(\alpha)$\;}

\Return $\Gamma^\alpha$.

\caption{Local conformal spatial prediction (LCSP).}\label{ag:lscp}
\end{algorithm}

With the help of the local exchangeability result in Theorem~\ref{thm:local.exchangeable}, we can establish the (asymptotically approximate) theoretical validity of the local spatial conformal prediction method.  As before, to highlight the conformal prediction interval's dependence on $s^\star$, we write it as $\Gamma^\alpha(Z^m; s^\star)$.

\begin{thm}
\label{thm:valid.local}
Consider an infill asymptotic regime with $n=N^2$ spatial grid points, with $N \to \infty$.  Fix an integer $m > 0$ and let the radius $r=r_N$ satisfy \eqref{eq:radius}.  Then under the assumptions of Theorem~\ref{thm:local.exchangeable}, the prediction intervals returned by  Algorithm LSCP are asymptotically valid at $s^\star$ in the sense that
\[ \inf \liminf_{N \to \infty} \prob^{m+1}\{ \Gamma^\alpha(Z^m; s^\star) \ni Y(s^\star) \} \geq 1-\alpha, \]
where $\prob^{m+1}$ is the joint distribution of $Y$ at the $m+1$ spatial grid points, as well as at $s^\star$, and the outermost infimum is over all such distributions $\prob$ for $Y$.  The inequality above can be made arbitrary close to equality by taking $m$ sufficiently large.
\end{thm}

\begin{proof}
See Appendix~\ref{proof:valid.local}.
\end{proof}

\subsection{Algorithm sLSCP: smoothed LSCP }

Theorem~\ref{thm:valid.local} implies that the local spatial conformal prediction with $m$ nearest neighbors is approximately valid under the infill asymptotic regime. However, in practice completely disregarding the contribution of the observations outside the $m$ nearest neighbors may be unsatisfactory, so we propose a smoothed version of the LSCP algorithm (sLSCP).

The GSCP algorithm weights all $n+1$ non-conformity measures $\delta_i$ equally in the plausibility computation in \eqref{eq:plausibility}, but this is questionable for non-stationary processes with  stochastic properties that vary throughout the spatial domain.  To allow for non-stationarity, the sLSCP algorithm weights the non-conformity measures $\delta_i$ by how far the corresponding $s_i$ is from the prediction location. Let $f$ be a non-increasing function, and define weights
\[ w_i \propto f(d_i), \quad i=1,\ldots,n+1, \]
where $d_i = \|s_i-s^\star\|$, $d_{n+1} \equiv 0$, and the proportionality constant is $\sum_{i=1}^{n+1} w_i$ so that $\sum_{i=1}^{n+1} w_i = 1$.  Different $f$ functions can be applied, but we recommend the normalized Gaussian kernel function with bandwidth $\eta$,
\begin{equation}
\label{eq:weight}
w_{i} = \frac{\exp(-d_{i}^2/2\eta^2)}{1+\sum_{j=1}^{n}\exp(- d_{j}^2/2\eta^2)}, \quad i=1,\ldots,n+1.
\end{equation}
Note that, if $\eta \to \infty$, then $w_i \to (n+1)^{-1}$ for each $i$, which corresponds to the GSCP algorithm.  Finally, with these new weights, the plausibility contour at a provisional value $y_{n+1}$ of $Y(s^\star)$ is given by
\begin{equation}\label{eq:sLSCP:plausibility}
  p(y_{n+1} \mid Z_{(n+1)}^{n+1}, w^{n+1}) = \sum_{i=1}^{n+1} w_i 1\{\delta_i \geq \delta_{n+1}\},
\end{equation}
As before, while many of non-conformity measures are possible, we recommend the Kriging-based strategy with $\delta_i=e_i^2$ and $e_i$ defined in \eqref{eq:residual}.

Since we are interested only in the local structure of $Y$, it is natural that locations far from $s^\star$ have negligible weight, as in \eqref{eq:weight}.  But including all $n$ observations requires some non-trivial calculations, e.g., inverting a large $n \times n$ covariance matrix.
Therefore, to avoid cumbersome and irrelevant computation, we recommend using only the $M << n$ observations closest to $s^\star$ for both the Kriging predictions that determine $\delta_i$ and in the plausibility scores in \eqref{eq:sLSCP:plausibility}.  The resulting method is both locally adaptive and computationally efficient for large data sets.

The tuning parameter $\eta$ can be selected using cross validation, as illustrated in Sections \ref{s:sim} and \ref{s:app}.  The value of $M$ is determined by the bandwidth $\eta$ so that all observations with substantial $w_i$ are included, as are observations that are required for the Kriging prediction of these observations. Typically the number of nearby observations required to approximate the Kriging prediction is a small subset of the total number of observations \citep{stein2002screening}.  As a rule of thumb, $M$ could be selected to roughly include all observations within $2\eta + r^\star$ radius of $s^\star$, where $2\eta$ captures observations with substantial weights, and $r^\star$ is selected so that all the $M$ observations include the nearest 15 neighbors of the observation within $2\eta$ of $s^\star$.
We summarize the details of Algorithm sLSCP in Algorithm~\ref{ag:lscp+}.
For simplicity, we use sLSCP and LSCP indistinguishably.

\begin{algorithm}[t]
\SetAlgoLined
\KwIn{observations $z_i = (s_i, y_i), i=1, \ldots, n$; predict location $s^\star$; non-conformity measure $\Delta$; significance level $\alpha$; and a fine grid of candidate $y^\star$ values}
\parameter{weight parameter $\eta \in (0, \infty)$; number of neighbors to consider $M \leq n$}
\KwOut{$(1-\alpha)100\%$ prediction interval, $\Gamma^\alpha$, for $Y(s^\star)$}

determine $M$ through $\eta$ if not given\;

form $z_i, i = 1, \ldots, M,$ based on $M$ locations closest to $s^\star$\;

$s_{M+1} \gets s^\star$\;

calculate weights $w_i, i=1, \ldots, M+1$ as in \eqref{eq:weight}\;

\For{$y_{M+1}$ in the specified grid}{

\For{i = 1 to M+1}{
define $z_{(i)}^{M+1}$ by removing $y_i$ from $z^{M+1}$\;
$\delta_i \gets \Delta(z_{(i)}^{M+1}, y_i)$ \;}

compute plausibility for $y_{M+1}$ as $p(y_{M+1} \mid z_{(M+1)}^{M+1}, w^{M+1})$ in \eqref{eq:sLSCP:plausibility}\;

include $y_{M+1}$ in $\Gamma^\alpha$ if $p(y_{M+1} \mid z_{(M+1)}^{M+1}, w^{M+1}) \geq t_M(\alpha)$\;}

\Return $\Gamma^\alpha$.

\caption{Smoothed local spatial conformal prediction (sLSCP).}\label{ag:lscp+}
\end{algorithm}

\section{Simulation study}\label{s:sim}

\subsection{Data generation} \label{s:sim:data}

We consider one mean-zero Gaussian stationary process (Scenario 1) and seven non-Gaussian and/or nonstationary data-generating scenarios (Scenarios 2 -- 8). The seven scenarios are generated based on transformations of a mean-zero stationary Gaussian process $X(s)$ and white noise process $E(s)$ with standard normal distribution, where $s=(s_x,s_y) \in [0,1]^2$. The stationary Gaussian process process $X(s)$ has a $\Matern$ covariance function with variance $\sigma^2=3$, range $\phi=0.1$ and smoothness $\kappa = 0.7$. Data are sampled on the $N\times N$ grid of $n=N^2$ points in the unit square, $s\in \{N^{-1}, 2 N^{-1}, \ldots, 1\}^2$, with $N = 20$ or $N = 40$.

Specifically, the eight scenarios under consideration are:
\begin{enumerate}
    \item \label{scen:stationary} $Y(s) = X(s) + E(s)$;
    \item \label{scen:cubic} $Y(s) = X(s)^3 + E(s)$;
    \item \label{scen:skewed} $Y(s) = q[\Phi\{X(s)/\sqrt{3}\}] + E(s)$ where $\Phi$ is the standard normal distribution function and $q$ is the $\text{Gamma}(1,3^{-1/2})$ quantile function;
    \item \label{scen:prod} $Y(s) = \sqrt{3} X(s)  |E(s)|$;
    \item \label{scen:gfun} $Y(s) = \mbox{sign}\{X(s)\}  |X(s)|^{s_x+1}+E(s)$;
    \item \label{scen:eastwest} $Y(s) = \sqrt{\omega(s) / 3}\, X(s) + \sqrt{1-\omega(s)} \, E(s)$ where $\omega(s)=\Phi(\frac{s_x-0.5}{0.1})$;
    \item \label{scen:nugget} $Y(s) = X(s) + s_x \, E(s)$; and
    \item \label{scen:spike} $Y(s) = X(s) + 10 \exp(-50\|s-c\|^2)$ where $c = (0.5,0.5)$.
\end{enumerate}
Scenario~\ref{scen:stationary} is Gaussian and stationary, Scenarios~\ref{scen:cubic}--\ref{scen:prod} are stationary but non-Gaussian, and Scenarios~\ref{scen:gfun}--\ref{scen:spike} are  nonstationary either in the spatial variance (Scenarios~\ref{scen:gfun} and \ref{scen:eastwest}), error variance (Scenario~\ref{scen:nugget}), or mean (Scenario~\ref{scen:spike}).  Figure~\ref{f:sim:overview} shows one realization of all the scenarios with $N=40$.

\begin{figure}\caption{One realization each from simulation scenarios with $N=40.$}\label{f:sim:overview}
\centering

\includegraphics[width=0.76\textwidth]{sim_scenarios_overview.pdf}
\end{figure}

\subsection{Prediction methods and metrics}\label{s:sim:methods}

For each dataset we apply the global and local (with $\eta=0.1$) conformal spatial prediction algorithms. For the parametric Kriging method and the initial Kriging predictions of our proposed conformal prediction, we estimate the spatial covariates parameters $\Theta$ (Section~\ref{s:back:spatial}) using empirical variograms methods \citep{cressie1992statistics}. The empirical variograms are calculated using the {\tt variog} function in the {\tt R} package {\tt geoR}, and the covariance parameters are chosen to minimize the weighted (by number of observations) squared error between the empirical and model-based variograms.

We compare the conformal methods with standard global Kriging prediction using \eqref{eq:kriging} to construct prediction interval, and the local Kriging (laGP) method of \cite{gramacy2015local} that dynamically defines the support of a Gaussian process predictor based on a local subset of the data. For laGP, we use the function provided by the {\tt laGP} package in {\tt R} the local sequential design scheme starting from 6 points to 50 points through an empirical Bayes' mean-square prediction error criterion.

Methods are compared using leave-one-out cross validation, i.e., we generate a prediction interval for each observation given all other observations. Each scenario is repeated 100 times, and performance is evaluated using average coverage of $(1-\alpha)100$\% prediction intervals, average interval width, and average interval score \citep{gneiting2007strictly}, defined as
$$S_\alpha(I_l,I_u;y^n) = \frac{1}{n}\sum_{i=1}^n\big[(I_u-I_l) + \tfrac{2}{\alpha} (I_l-y_i) 1(y_i<I_l) + \tfrac{2}{\alpha}(y_i-I_u)1(y_i>I_u)\big],$$ where $[I_l,I_u]$ represents the $100(1-\alpha)\%$ prediction interval and $y^n=(y_1,\ldots,y_n)$ are the observations.  A smaller interval score is desirable as this rewards both high coverage and narrow intervals. We use $\alpha=0.1$ in this simulation study.

\subsection{Results}\label{s:sim:results}

We present the results for Scenarios~\ref{scen:stationary}--\ref{scen:eastwest} in this section; the others are in Appendix \ref{s:A3}. For the stationary Scenarios~\ref{scen:stationary}--\ref{scen:prod} we present results averaged over data sets and all spatial locations in Table~\ref{t:sim:stationary}. For the non-stationary scenarios varying across $s_x$ (Scenarios \ref{scen:gfun}--\ref{scen:nugget}), we present the results by the first spatial coordinate ($s_x$) averaged over the data sets and the second coordinate ($s_y$) in Figures~\ref{f:gfun_n_eastwest}, e.g., the value of coverage plotted at $s_x=N^{-1}$ is the average of the coverage over the $N$ points of the form $(N^{-1},s_y)$ for $s_y\in\{N^{-1}, 2 N^{-1}, \ldots, 1\}$.

\begin{table}[t]
\caption{Performance comparison for simulation Scenarios (``Scen'') \ref{scen:stationary} -- \ref{scen:prod}. The metrics are the empirical coverage of 90\% prediction intervals (``Cov90''), the width of prediction intervals (``Width'') and the interval score (``IntScore''), each averaged over location and dataset. The methods are global (GSCP) and local (LSCP) conformal prediction, stationary and Gaussian Kriging (``Kriging'') and local approximate Gaussian process (``laGP'') regression. }\vspace{6pt} \label{t:sim:stationary}
\centering
\begin{tabular}{ll|ccc|ccc}
 & & \multicolumn{3}{c|}{$N=20$} & \multicolumn{3}{c}{$N=40$}   \\
 Scen & Method & Cov90 & Width & IntScore & Cov90 & Width & IntScore \\
  \hline
\ref{scen:stationary}
  & GSCP    & 89.9\% & 4.58 & 5.77 & 90.0\% & 3.99 & 5.01 \\
  & LCSP    & 88.5\% & 4.55 & 5.96 & 89.2\% & 3.94 & 5.07 \\
  & Kriging	& 90.7\% & 4.68 & 5.75 & 88.5\% & 3.82 & 5.01 \\
  & laGP	& 86.9\% & 4.74 & 6.52 & 88.7\% & 4.26 & 5.52 \\
&&&&&&&\vspace{-6pt}\\
\ref{scen:cubic}
  & GSCP    & 89.9\% & 31.84 & 63.07 & 90.0\% & 21.65 & 41.43 \\
  & LCSP    & 90.7\% & 30.72 & 55.32 & 90.9\% & 21.06 & 33.77 \\
  & Kriging	& 93.7\% & 43.70 & 65.54 & 92.4\% & 26.77 & 42.47 \\
  & laGP	& 91.7\% & 40.07 & 67.49 & 93.4\% & 32.07 & 46.27 \\
&&&&&&&\vspace{-6pt}\\
\ref{scen:skewed}
  & GSCP    & 89.9\% & 4.64 & 6.22 & 90.0\% & 4.00 & 5.14 \\
  & LCSP    & 88.7\% & 4.58 & 6.22 & 89.2\% & 3.94 & 5.14 \\
  & Kriging	& 90.9\% & 4.78 & 6.19 & 88.9\% & 3.87 & 5.15 \\
  & laGP	& 87.8\% & 4.73 & 6.71 & 89.0\% & 4.27 & 5.61 \\
&&&&&&&\vspace{-6pt}\\
\ref{scen:prod}
  & GSCP    & 89.9\% & 6.93 & 11.03 & 90.0\% & 6.21 & 9.84 \\
  & LCSP    & 88.9\% & 6.82 & 11.21 & 90.0\% & 5.96 & 9.26 \\
  & Kriging	& 91.7\% & 7.65 & 11.02 & 91.2\% & 6.61 & 9.85 \\
  & laGP	& 89.2\% & 7.37 & 11.93 & 91.0\% & 6.68 & 10.08 \\
\end{tabular}
\end{table}

Table~\ref{t:sim:stationary} shows the results for data generated from stationary processes (Scenarios \ref{scen:stationary}--\ref{scen:prod}).

In Scenario~\ref{scen:stationary}, the Gaussian and stationary process, the performance of GSCP, LSCP, and Kriging are comparable. Kriging performs well in this case since the data generating mechanism aligns with its underlying assumption, but the conformal methods are competitive with the parametric model.

In Scenarios~\ref{scen:cubic}, \ref{scen:skewed}, and \ref{scen:prod}, the non-Gaussian but stationary processes, GSCP, LSCP, and Kriging perform more or less the same in terms of interval score and outperform laGP. However, the coverage of the conformal methods, especially the GSCP algorithm, is closer to the nominal level than Kriging and the Kriging intervals are generally wider than the conformal intervals.  Figure~\ref{f:skew_variance} examines the results in more detail for Scenario 3.  Since this is a stationary process, the width of the prediction intervals should be similar across different locations, which is the case for the global methods GSCP and Kriging but not the local methods LSCP and laGP. This extra variation is reflected in the coverage proportion across datasets. The coverage of GSCP stabilizes at 90\% for each simulation while the coverage for other methods varies considerably and is often below the desired level.  This detailed analysis focused on Scenario~3, but the same conclusions can be drawn in Scenarios~\ref{scen:cubic} and \ref{scen:prod}.

\begin{figure}
\centering

\includegraphics[width=1\textwidth]{sim_skew_variance.pdf}
\caption{Detailed performance comparison for Scenario 3 with $Y(s) = q[\Phi\{X(s)/\sqrt{3}\}] + E(s)$.  For each dataset and method we compute the standard deviation of interval width across space and the average coverage across space, and the boxplots summarize the sampling distribution of these statistics across the 100 simulated datasets.}
\label{f:skew_variance}
\end{figure}

Figure \ref{f:gfun_n_eastwest} shows the results for nonstationary Scenarios \ref{scen:gfun} and \ref{scen:eastwest} when $N=40$. The results for $N=20$ and the results of the other nonstationary scenarios are in Appendix \ref{s:A3}.  LSCP performs the best among the four methods for these nonstationary scenarios. For Scenario \ref{scen:gfun}, the process is Gaussian to the west and non-Gaussian to the east. The global prediction methods GSCP and Kriging generate prediction intervals with similar width for all $s_x$ (ignoring edge effects), while LSCP and laGP provide wider intervals on the east ($s_x$ near 1) than the west ($s_x$ near 0). LSCP has coverage around 90\% for $s_x$, and the lowest interval scores, especially in the east where the process is more non-Gaussian.  Similarly, in Scenario \ref{scen:eastwest}, the correlation is stronger in the east than the west, and the LSCP performs the best by providing adaptive prediction interval width and valid coverage across space.

\begin{figure}
\centering

\includegraphics[width=0.9\textwidth]{sim_result_gfun_n_eastwest_40.pdf}
\caption{Performance comparison for Scenario \ref{scen:gfun}: $Y(s) = \mbox{sign}\{X(s)\} \cdot |X(s)|^{s_x+1}+E(s)$ and Scenario \ref{scen:eastwest}: $Y(s) = \sqrt{\omega(s)}\cdot X(s)/\sqrt{3} + \sqrt{1-\omega(s)}\cdot E(s)$ where $\omega(s)=\Phi(\frac{s_x-0.5}{0.1})$ when $N=40$. }
\label{f:gfun_n_eastwest}
\end{figure}

We also conducted a simulation study when observations are sampled uniformly on $[0,1]^2$. The performance is very similar to that sampled on equally-spaced grid points so we only show the results of sampling on grid points in this paper.

\subsection{Sensitivity Analysis}\label{s:sim:sensitivity}

Finally, we study sensitivity of the proposed method to tuning parameter selection and to the initial estimates of the $\Matern$ covariance parameters. Figure \ref{f:sensitivity} gives the average interval score of LSCP with different tuning parameter $\eta$ for Scenarios \ref{scen:gfun} and \ref{scen:eastwest}, along with the average interval score of  GSCP, Kriging, and laGP. We find that there is a V-shaped relationship between resulting interval score and tuning parameter $\eta$. The average performance is the best when $\eta=0.2$, and LSCP outperforms the other methods regardless the selection of $\eta$. To achieve the best performance, in practice, we would select $\eta$ using cross validation.

\begin{figure}
\centering

\includegraphics[width=1\textwidth]{sim_sensitivity.pdf}
\caption{Sensitivity analysis for simulation Scenarios \ref{scen:gfun} and \ref{scen:eastwest}.  The interval score for the LSCP method is computed for several kernel bandwidth parameters, $\eta$, and compared to other methods.}
\label{f:sensitivity}
\end{figure}

The proposed spatial conformal prediction methods rely on estimation of the $\Matern$ covariance parameter, $\Theta$. To test for sensitivity to this step, we repeat simulation Scenario 1 since in this case the optimal $\Theta$ is known to be the true $\Matern$ covariance parameters used to generate the data. When generating prediction intervals, rather than estimating the $\Matern$ covariance parameter $\hat\Theta$ from the simulated data, we use the values of $\hat\Theta$ listed in the first four columns of Table \ref{t:sensitivity_thetaHat} for conformal prediction and Kriging, where each of the four parameters is modified by plus or minus 50\% of the true value. We use the GSCP method, $\alpha=0.1$ and $N=20$ for this sensitivity analysis. The corresponding average coverage and average width of the conformal prediction intervals and the Kriging intervals are listed in the last four columns.

\begin{table}
\centering
\caption{Sensitivity analysis for the simulation study. The first four columns represent the $\Matern$ covariance parameter ($\Theta$) used to generate the conformal and Kriging prediction intervals; the first row gives the true value. The last four columns give the average coverage and interval width of the 90\% intervals produced by the conformal and Kriging methods when using the corresponding $\Theta$.} \label{t:sensitivity_thetaHat}\vspace{6pt}
\begin{tabular}{cccc|cc|cc}
Nugget& Partial sill & Range  & Smoothness & \multicolumn{2}{c|}{GSCP} & \multicolumn{2}{c}{Kriging} \\
($\tau^2$) &($\sigma^2$) &($\phi$)& ($\kappa$) & Cov90 & Width & Cov90 & Width \\
  \hline
1.0 & 3.0 & 0.10 & 0.70 & 89.69\% & 4.57 & 90.05\% & 4.60 \\
1.5 & 3.0 & 0.10 & 0.70 & 89.69\% & 4.61 & 93.73\% & 5.24 \\
0.5 & 3.0 & 0.10 & 0.70 & 89.70\% & 4.59 & 82.76\% & 3.84 \\
1.0 & 4.5 & 0.10 & 0.70 & 89.68\% & 4.63 & 92.64\% & 5.04 \\
1.0 & 1.5 & 0.10 & 0.70 & 89.68\% & 4.63 & 85.25\% & 4.09 \\
1.0 & 3.0 & 0.15 & 0.70 & 89.68\% & 4.60 & 86.71\% & 4.21 \\
1.0 & 3.0 & 0.05 & 0.70 & 89.67\% & 4.61 & 95.13\% & 5.54 \\
1.0 & 3.0 & 0.10 & 1.05 & 89.69\% & 4.58 & 85.56\% & 4.08 \\
1.0 & 3.0 & 0.10 & 0.35 & 89.70\% & 4.63 & 94.98\% & 5.53 \\
\end{tabular}
\end{table}

While using the true $\Theta$ leads to slightly more narrow intervals, the average coverage of the conformal prediction intervals are between 89.67\% and 89.70\% across $\Theta$.  In contrast, the coverage of the Kriging intervals ranges from 82.76\% to 95.13\% across $\Theta$.  Therefore, the conformal methods are more robust to estimation of $\Matern$ covariance parameters than the parametric Kriging model.

\section{Real data analysis}\label{s:app}

This section demonstrates the performance of conformal prediction method using the canopy height data in Figure~\ref{f:canopy_orig}. The data were originally presented in \cite{cook2013nasa} and were analyzed using a stationary Kriging model in \cite{datta2016hierarchical}. The data are available in the {\tt spNNGP} package \citep{finley2017spnngp} in {\tt R}. There are $n=1,723,137$ observations and clear nonstationarity and non-normality.  For example, there are several heterogeneous areas with small height canopies around the location with longitude and latitude being 729,000 and 470,000, respectively.

We compare methods using 90\% prediction intervals for 10,000 randomly chosen test locations. We use LSCP and choose the tuning parameters as $\eta=6\times10^{-4}$ because this minimized the interval score averaged over 1,000 randomly chosen validation locations. The average interval score and tuning parameter $\eta$ also show a V-shape relationship similar to that in Figure~\ref{f:sensitivity}. Table~\ref{tab:canopy_perf_comp} compares the performance of LSCP on the 10,000 test locations with Kriging and laGP.  LSCP outperforms the other methods as the empirical coverage of LSCP is the closet to the desired 90\% and the LSCP minimizes the interval score.

\begin{table}
    \centering
    \caption{Performance comparison for the canopy height data. The metrics are the width, coverage (``Cov90'') and interval score (``IntScore'') of 90\% prediction intervals, each averaged over 10,000 randomly chosen test locations. The methods are local  conformal prediction (LSCP), stationary and Gaussian Kriging (``Kriging'') and local approximate Gaussian process (``laGP'') regression.}\vspace{6pt}
    \begin{tabular}{l|ccc}
             & Width & Cov90   & IntScore\\
     \hline
     LSCP    & 3.00  & 89.2\%  & 4.31 \\
     Kriging & 5.54  & 96.7\%  & 6.89 \\
     laGP    & 5.04  & 91.3\%  & 6.61 \\
    \end{tabular}
    \label{tab:canopy_perf_comp}
\end{table}

Figure \ref{f:canopy} plots the interval widths for each method. Unlike Kriging, the LSCP and laGP interval widths are locally adaptive with wider intervals in heterogeneous areas and more narrow intervals in homogeneous areas. Comparing LSCP and laGP, LSCP generally provides narrower intervals than laGP, which means the proposed method is more efficient than laGP. In addition, the locations of the observations that fall outside the prediction intervals are uniformly distributed for LSCP and laGP, but clustered in the heterogeneous areas for Kriging.  In short, the proposed spatial conformal prediction algorithm shows its superiority in this real data analysis.

\begin{figure}
  \begin{subfigure}[b]{0.5\linewidth}
    \centering

\includegraphics[width=0.9\linewidth]{canopy_result_summary_part1.pdf}
    \caption{Original data}
    \label{f:canopy_orig}
  \end{subfigure}
  \begin{subfigure}[b]{0.5\linewidth}
    \centering

\includegraphics[width=0.9\linewidth]{canopy_result_summary_part2.pdf}
    \caption{LSCP result}
    \label{f:canopy_LSCP}
  \end{subfigure}
  \newline
  \begin{subfigure}[b]{0.5\linewidth}
    \centering

\includegraphics[width=0.9\linewidth]{canopy_result_summary_part3.pdf}
    \caption{Kriging result}
    \label{f:canopy_Kriging}
  \end{subfigure}
  \begin{subfigure}[b]{0.5\linewidth}
    \centering

\includegraphics[width=0.9\linewidth]{canopy_result_summary_part4.pdf}
    \caption{laGP result}
    \label{f:canopy_laGP}
  \end{subfigure}
\caption{(a) Heatmap of the canopy height data; (b)--(d) Prediction interval width (color) and locations not covered (points) for LSCP, Kriging, and laGP. Longitude and latitude are in UTM Zone 18. }
\label{f:canopy}
\end{figure}

\section{Discussion}\label{s:discussion}

In this paper, we proposed a spatial conformal prediction algorithm to provide valid, robust, and model-free prediction intervals by combining spatial methods and the classical conformal prediction. We provided both global and local versions to accommodate different stationarity cases and sampling designs. We proved their validity under various sampling designs and data-generating mechanisms. To the authors' knowledge, this work is the first in making the classical conformal algorithm work for non-exchangeable data. Our simulation studies and real data analyses demonstrate the advantage of the proposed spatial conformal prediction algorithms. We also developed an {\tt R} package entitled {\tt scp} (\url{https://github.com/mhuiying/scp}) to compute plausibility and generate spatial prediction intervals using either Kriging residual or any other user-defined non-conformity measure.

An attractive feature of the proposed algorithms is that they are model-free in the sense that their theoretical validity does not depend on specifying the correct model. In our implementations we use the squared residuals from a simple parametric model to define the non-conformity terms. However, as anticipated by our theoretical results, our simulation study shows that the methods work well even if the parametric family or the mean and covariance functions are misspecified, and that the results are insensitive to inaccurate estimation of the parameters in the parametric model. This robustness allows the methods to be applied broadly and with confidence.

The local conformal prediction method relies on a dense grid of points around each prediction location and thus a large dataset.  Computation time often prohibits application of spatial methods to large datasets.  Fortunately, we are able to apply the conformal prediction method to large datasets by exploiting local algorithms and an explicit solution of the plausibility contour assuming the non-conformity measure is the squared Kriging residual. Similar derivations are needed for prediction procedures other than Kriging in order to maintain computational efficiency.

Future research directions include extending the work to spatial processes with discrete observations, e.g., when $Y(s)$ is binary or a count.  Generalized spatial linear models are compatible with our current framework, see \eqref{eq:glm}, but continuity in the distribution function is required.  Therefore, further studies would be required to establish the validity of conformal prediction for discrete data.  It is also of interest to extend the proposed spatial conformal prediction methods to temporal or multivariate spatial data.

\section*{Acknowledgments}

HM (DMS--1638521) and RM (DMS--1811802) are supported by the National Science Foundation. BJR is supported by the National Institutes of Health (R01ES031651 and R01ES027892) and King Abdullah University of Science and Technology (3800.2).

\begin{singlespace}
	\bibliographystyle{apalike}
	\bibliography{refs}
\end{singlespace}

\appendix

\section{Computational details}\label{s:A1}

Applying the predictive distributions in \eqref{eq:kriging} and provisionally setting $Y_{n+1}=y$, for $i\ne n+1$ we can write
\begin{equation}\label{e:quad}
\delta_i-\delta_{n+1} = q_{ii}\left(Y_i-{\tilde Y}_{i,n+1}+\frac{q_{i,n+1}}{q_{ii}}y\right)^2 - q_{n+1,n+1}\left(y-{\hat Y}_{n+1}\right)^2 = U_{i} + V_{i}y+W_{i}y^2
\end{equation}
where ${\tilde Y}_{i,n+1} = -\sum_{j\ne\{i,n+1\}}q_{ij}Y_j/q_{ii}$, $U_{i} = q_{ii}(Y_i-{\tilde Y}_{i,n+1})^2-q_{n+1,n+1}{\hat Y}_n^2$, $V_{i} = 2q_{i,n+1}(Y_i-{\tilde Y}_{i,n+1}) + 2q_{n+1,n+1}{\hat Y}_{n+1}$ and $W_{i} = q_{i,n+1}^2/q_{ii} -q_{n+1,n+1}$.

Since the inverse covariance matrix $Q$ is positive definite, $q_{ii}>0$ and $q_{ii}>\max_{j\ne i}|q_{ij}|$, so $W_{i} <0$ and $V_i^2-4U_iW_i = \frac{4q_{n+1,n+1}}{q_{ii}}[\sum_{i=1}^n q_{ij}Y_j+q_{i,n+1}\hat Y_{n+1}]^2\geq0$. Therefore, the $y$ satisfying $\delta_i-\delta_{n+1}\geq0$ are within two roots, denoted as $a_i \leq b_i$, of the quadratic equation $U_{i} + V_{i}y+W_{i}y^2=0$. Then the plausibility calculation \eqref{eq:GSCP:plausibility} becomes $p(y|Z_{(n+1)}^{n+1})=\frac{1}{n+1}\sum_{i=1}^{n}1\{U_{i} + V_{i}y+W_{i}y^2\geq0\}+\frac{1}{n+1} = \frac{1}{n+1}\sum_{i=1}^{n}1\{a_i \leq y \leq b_i\}+\frac{1}{n+1}$, which is a step function with jumping points being $a_i$'s and $b_i$'s, and the plausibility function for sLSCP in \eqref{eq:sLSCP:plausibility} simply becomes
$\sum_{i=1}^{n}w_i1\{a_i \leq y \leq b_i\}$.

With the help of the step function, we can solve for the prediction interval directly with no need to enumerate for possible $y$ satisfying $p(y|Z_{(n+1)}^{n+1})\geq t_n(\alpha)$.

\section{Proofs}

\subsection{Proof of Proposition~\ref{prop:sphere}}
\label{proof:sphere}

To be even more explicit about notation, for a given set of data $\{s^n, y^n\}$ and a new location $s_{n+1}$, write $\Gamma^\alpha(s^n, y^n; s_{n+1})$.  Now define the conditional coverage probability function
\[ c_\alpha(s^\star) = \prob^{n+1}\{\Gamma^\alpha(S^n, Y^n; S_{n+1}) \ni Y_{n+1} \mid S_{n+1} = s^\star\}, \]
where $\prob^{n+1}$ is the joint distribution of $(S_1,Y_1),\ldots,(S_{n+1}, Y_{n+1})$.  Let $\rho$ denote a generic $3 \times 3$ rotation matrix.  From our Kriging-based construction of the conformal prediction interval, it is each to check that $\Gamma^\alpha(\rho s^n, y^n; \rho s_{n+1}) = \Gamma^\alpha(s^n, y^n; s_{n+1})$ for all $\rho$ and all $(s^n, y^n, s_{n+1})$, where $\rho s^n$ denotes the rotation applied to each location in $s^n$.  Moreover, by the uniformity and isotropic stationarity assumptions, respectively, the marginal distribution of $S^{n+1}$ and the conditional distribution of $Y^{n+1}$, given $S^{n+1}$, are invariant to rotations of the locations in the sphere.  Therefore, $c_\alpha(\rho s^\star) \equiv c_\alpha(s^\star)$, i.e., the conditional coverage probability is constant in $s^\star$.  But the result in Theorem~\ref{thm:valid.global} applies in this scenario as well, so if the marginal coverage probability exceeds $1-\alpha$, then the constant conditional coverage probability must exceed $1-\alpha$ too.

\subsection{Proof of Theorem~\ref{thm:local.exchangeable}}
\label{proof:local.exchangeable}

Write the localized process as
\[ \tilde Y_r(u) = \psi \bigl( \tilde X_r(u), \tilde E_r(u) \bigr), \quad u \in \calU, \]
where $\psi$ is a continuous real-valued function of two arguments.  We are interested in the finite-dimensional distribution of the localized process, so take a fixed set of $m$ vectors $u_1,\ldots,u_m$ in $\calU$.  Define the vectors
\[ \tilde X_r^m = \bigl( \tilde X_r(u_1),\ldots, \tilde X_r(u_m) \bigr) \quad \text{and} \quad \tilde E_r^m = \bigl( \tilde E_r(u_1), \ldots, \tilde E_r(u_m) \bigr). \]
First, by Markov's inequality, we have that
\[ \prob\{\|\tilde X_r^m - X(s^\star) \, 1_m\| > \delta\} \leq \delta^{-2} \, \E\|\tilde X_r^M - X(s^\star)\, 1_m\|^2, \quad \text{for any $\delta > 0$}, \]
where $\|\cdot\|$ denotes the usual norm on $\RR^m$ and $1_m$ is a $m$-vector of unity.  The expectation in the upper bound can be rewritten as
\begin{align*}
\E\|\tilde X_r^m - X(s^\star) \, 1_m\|^2 & = \sum_{i=1}^m \E\{\tilde X_r(u_i) - X(s^\star)\}^2 \\
& = \sum_{i=1}^m \E\{X(s^\star + r u_i) - X(s^\star)\}^2,
\end{align*}
and, since $m$ is fixed, the assumed $L_2$-continuity of the $X$ process implies that the right-hand side vanishes as $r \to 0$.  Therefore, we have that $\tilde X_r^m \to X(s^\star) 1_m$ in probability and, hence, in distribution.  For the $E$ part in the above representation of $\tilde Y_r$, the locally iid assumption implies
\[ \tilde E_r^m \to \tilde E_0^m, \quad \text{in distribution, as $r \to 0$} \]
where $\tilde E_0^m$ is an iid vector.  This, together with the continuous mapping theorem, gives
\[ \tilde Y_r^m \to \psi\bigl( X(s^\star) \, 1_m, \tilde E_0^m \bigr) \quad \text{in distribution, as $r \to 0$}, \]
where $\psi$ is being applied component-wise.  The right-hand side is clearly conditionally iid, given $X(s^\star)$, which implies exchangeability.  Finally, since this holds for all $m$ and for all $(u_1,\ldots,u_m)$, we know that the process $\tilde Y_r$ has a distributional limit, and that limit is an exchangeable process.

\subsection{Proof of Theorem~\ref{thm:valid.local}}
\label{proof:valid.local}

For the $m$ spatial grid locations $s_{i_j} = s^\star + r u_j$, for $j=1,\ldots,m$, in the radius-$r$ neighborhood around $s^\star$, write $\tilde Y_r^{m+1}$ for the $(m+1)$-vector that includes $Y(s^\star + r u_j)$, for $j=1,\ldots,m$, and $Y(s^\star)$ appended as the $(m+1)^\text{st}$ value.  Let $\delta_j$, for $j=1,\ldots,m,m+1$, denote the corresponding non-conformity scores.  Collect these together into the $(m+1)$-vector $\delta_r^{m+1}$.  Then
\[ Y(s^\star) \in \Gamma^\alpha(Z^m; s^\star) \iff \rank(m+1; \delta_r^{m+1}) > \lfloor (m+1) \alpha \rfloor, \]
where $\rank(m+1; \delta_r^{m+1})$ denotes the rank of the $(m+1)^\text{st}$ entry within the collection $\delta_r^{m+1}$.  By Theorem~\ref{thm:local.exchangeable}, we have that $\tilde Y_r^{m+1}$ converges in distribution, as $r \to 0$ or, equivalently, as $N \to \infty$, to an exchangeable $\tilde Y_0^{m+1}$.  And since our Kriging-based non-conformity measure is continuous, so that $\delta_r^{m+1}$ has a density with respect to Lebesgue measure on $\RR^{m+1}$, it follows easily from the Portmanteau lemma \citep[e.g.,][Lemma~2.2]{vaart1998}, that $\delta_r^{m+1}$ converges in distribution to $\delta_0^{m+1}$, say, which is just the non-conformity measure applied to the limit $\tilde Y_0^{m+1}$.  Similarly, it is easy to show that
\[ \rank(m+1; \delta_r^{m+1}) \to \rank(m+1; \delta_0^{m+1}) \quad \text{in distribution, as $r \to 0$ or $N \to \infty$}. \]
Since the limit $\tilde Y_0^{m+1}$ is exchangeable, the structure of the non-conformity measure, implies that $\delta_0^{m+1}$ is also exchangeable.  Therefore, $\rank(m+1; \delta_r^{m+1})$ converges in distribution to a random variable $U$, uniformly distributed on $\{1,\ldots,m,m+1\}$, so
\[ \prob^{m+1}\{ Y(s^\star) \in \Gamma^\alpha(Z^m; s^\star) \} \to \prob\{ U > \lfloor (m+1)\alpha \rfloor\} \geq 1-\alpha, \quad \text{as $r \to 0$ or $N \to \infty$}. \]
Finally, the probability on the right-hand side of ``$\to$'' converges from above to $1-\alpha$ as $m \to \infty$, so the inequality can be made arbitrarily close to equality with $m$ sufficiently large.

\section{Additional simulation results}
\label{s:A3}

This appendix presents the additional simulation results for Section \ref{s:sim}. Figure \ref{f:gfun_n_eastwest_20} shows the performance of GSCP, LSCP, Kriging, and laGP of the nonstationary Scenarios \ref{scen:gfun} and \ref{scen:eastwest} when $N=20$. The patterns are similar to what stated in the main text when $N=40$, and so are the conclusions. LSCP performs the best by giving the smallest interval scores, the desired coverage around 90\%, and adaptive interval width. LSCP is also the winning method for nonstationary Scenarios \ref{scen:nugget} and \ref{scen:spike}, and Figure \ref{f:nugget} ($N=20, 40$) and Figure \ref{f:spike} ($N=20, 40, 60$) show their results, respectively. For Scenario \ref{scen:spike}, the spike is injected at the center of $[0,1]^2, (0.5,0.5)$. The process pattern varies according to the distance to the center, so we present results by averaging over datasets and locations with the same distance to center.

\begin{figure}
\centering

\includegraphics[width=0.9\textwidth]{sim_result_gfun_n_eastwest_20.pdf}
\caption{Performance comparison for Scenario \ref{scen:gfun}: $Y(s) = \mbox{sign}\{X(s)\} \cdot |X(s)|^{s_x+1}+E(s)$ and Scenario \ref{scen:eastwest}: $Y(s) = \sqrt{\omega(s)}\cdot X(s)/\sqrt{3} + \sqrt{1-\omega(s)}\cdot E(s)$ where $\omega(s)=\Phi(\frac{s_x-0.5}{0.1})$ when $N=20$. }
\label{f:gfun_n_eastwest_20}
\end{figure}

\begin{figure}
\centering

\includegraphics[width=0.9\textwidth]{sim_result_nugget.pdf}
\caption{Performance comparison for Scenario \ref{scen:nugget}: $Y(s) = X(s) + \mbox{Normal}(0,s_x)$ when $N=20$ and $N=40$.}
\label{f:nugget}
\end{figure}

\begin{figure}
\centering

\includegraphics[width=0.9\textwidth]{sim_result_spike.pdf}
\caption{Performance comparison for Scenario \ref{scen:spike}: $Y(s) = X(s) + 10\cdot \exp(50||s-c|||^2)$ when $N=20, 40$ and 60.}
\label{f:spike}
\end{figure}

\end{document}

config: arxiv


TITLE: stuck on deriving Cauchy product
QUESTION [0 upvotes]: Could someone give me some heads up help or hint as to how I go from the LHS of the inequality to the other?$$\sum\limits_{k = n + 1}^{2n} {\sum\limits_{j = k - n}^n {{c^2}{r^{ - k}}{{\left| {z - {z_0}} \right|}^k} \le } } \sum\limits_{k = n + 1}^{2n} {(2n - k + 1){c^2}{q^k} \le {n^2}{c^2}{q^{n + 1}}}$$ where $q: = \left| {z - {z_0}} \right|/r \prec 1$. I believe that I understand how the part in the parentheses, () and the $q^k$ come about but don't get the last inequality.
$$\begin{array}{l}\text{I see the sum goes from }k - n\text{ to }n\text{ so}\\\text{there are }n - k + n + 1\text{ terms, namely }\\\text{there are  }2n - k + 1\text{ terms which is in the parentheses}\text{.}\\\text{The outer sum working on this}\end{array}$$
$$\begin{array}{l}\sum\limits_{n + 1}^{2n} {} \text{ there are n terms}\text{. }\\\sum\limits_{k = n + 1}^{2n} {(2n - k + 1)}  = 2\sum\limits_{}^{} n  - \sum\limits_{}^{} k  + \sum\limits_{}^{} 1  = 2{n^2} - n - \sum\limits_{}^{} k \\ \le 2{n^2}\end{array}$$
It comes from here, Proof of Cauchy product: Elias Wegert (auth.) - Visual Complex Functions_ An Introduction with Phase Portraits (2012, Birkhäuser Basel)

REPLY [0 votes]: $\begin{array}{l}\text{Since }\sum\limits_{k = n + 1}^{2n} {} \text{ : implies }k \ge n + 1\\\text{Also since }q \prec 1\text{ (hyp)}\\{q^k} \le {q^{n + 1}}\\\sum\limits_{k = n + 1}^{2n} {} \text{ has n terms (i}\text{.e}\text{. 2n - n - 1 + 1)}\\\text{therefore, }\sum\limits_{k = n + 1}^{2n} {(2n - k + 1)}  \le n.n = {n^2}\\\text{since, }\\(k \ge n + 1 \Leftrightarrow  - k \le  - n - 1 \Leftrightarrow 2n - k + 1 \le 2n - n = n)\end{array}
%
$


TITLE: For each $X\in G_{e}$, exists a unique smooth homomorphism $\phi:\mathbb{R}\to G$ with $\phi'(0)=X$ (Spivak)
QUESTION [2 upvotes]: Hello. I don't quite understand the proof of Corollary 8 in a comprehensive introduction to differential geometry by Spivak. If $\Phi:\mathbb{R}\to G_{e}=\mathcal{L}(G)$, by Theorem 5, exists an smooth homomorphism $\phi:U\to H$ such that ${\phi_{*}}_e\left(\left.\frac{d}{dt}\right|_{t=0}\right)=X$ but shouldn't ${\phi_{*}}_e\left(\left.\frac{d}{dt}\right|_{t=0}\right)=\Phi\left(\left.\frac{d}{dt}\right|_{t=0}\right)$? although this does not make sense because it should be a real number because the definition of $\Phi:\mathbb{R}\to G_e$ What am I misunderstanding?

REPLY [0 votes]: I just solved it. Is being identified $\alpha\in\mathbb{R}$ with $\alpha\left.\frac{d}{dt}\right|_{t=0}$ ($\alpha$ is the coordinate in the basis $\left.\frac{d}{dt}\right|_{t=0}$ of $T_{t=0}(\mathbb{R})$).


text:
\section{Numerical experiments}
\label{sec:experiment}
In this section, we validate our theory and demonstrate the efficiency of \fedac via experiments.\footnote{Code repository link: \url{https://github.com/hongliny/FedAc-NeurIPS20}.}
The performance of \fedac is tested against \fedavg (a.k.a., Local SGD), (distributed) Minibatch-SGD (\mbsgd) and Minibatch-Accelerated-SGD (\mbacsgd) \citep{Dekel.Gilad-Bachrach.ea-JMLR12,Cotter.Shamir.ea-NeurIPS11} on $\ell_2$-regularized logistic regression for UCI \texttt{a9a} dataset \citep{Dua.Graff-17} from \texttt{LibSVM} \citep{Chang.Lin-11}.
The regularization strength is set as $10^{-3}$.
The hyperparameters $(\gamma, \alpha, \beta)$ of \fedac follows \fedaci where strong-convexity $\mu$ is chosen as regularization strength $10^{-3}$.
We test the settings of $M = 2^2, \ldots, 2^{13}$ workers and $K = 2^0, \ldots, 2^8$ synchronization interval.
For all four algorithms, we tune the learning-rate $\eta$ \emph{only} from the same set of levels within $[10^{-3}, 10]$.
We choose $\eta$ based on the best suboptimality.
We claim that the best $\eta$ lies in the range $[10^{-3}, 10]$ for all algorithms under all settings.
We defer the rest of setup details to \cref{additional:expr}.
In \cref{fig:a9a:1e-3:M}, we compare the algorithms by measuring the effect of linear speedup
under variant $K$.
\begin{figure}[ht]
    \centering
    \includegraphics[width=\textwidth]{fig/a9a_1e-03_M_6blocks.pdf}
    \caption{\textbf{Observed linear speedup with respect to the number of workers $M$ under various synchronization intervals $K$.}
    Our \fedac is tested against three baselines \fedavg, \mbsgd, and \mbacsgd.
    While all four algorithms attain linear speedup for the fully synchronized ($K=1$) setting, \fedavg and \mbsgd lose linear speedup for $K$ as low as 8.
    \mbacsgd is comparably better than the other two baselines but still deteriorates significantly for $K \geq 64$.
    \fedac is most robust to infrequent synchronization and outperforms the baselines by a margin for $K \geq 64$.
    }
    \label{fig:a9a:1e-3:M}
\end{figure}

In the next experiments, we provide an empirical example to show that the direct parallelization of standard accelerated SGD may indeed suffer from instability.
This complements our \cref{instability:sketch}) on the initial-value instability of standard AGD.
Recall that \fedaci \cref{fedaci} and \fedacii \cref{fedacii} adopt an acceleration-stability tradeoff technique that takes $\gamma = \max \{ \sqrt{\frac{\eta}{\mu K}} , \eta\}$. Formally, we denote the following direct acceleration of \fedac without such tradeoff as ``vanilla \fedac'':
$
    \eta \in (0, \frac{1}{L} ], \gamma =  \sqrt{\frac{\eta}{\mu}} , \alpha  = \frac{1}{\gamma \mu},   \beta = \alpha + 1.
$
In \cref{fig:instability}, we compare the vanilla \fedac with the (stable) \fedaci and the baseline \mbacsgd.
\begin{figure}[ht]
    \centering
    \includegraphics[width=\textwidth]{fig/a9a_1e-04_M_6blocks.pdf}
    \caption{
        \textbf{Vanilla \fedac versus (stable) \fedaci and baseline \mbacsgd on the observed linear speedup w.r.t. $M$ under various synchronization intervals $K$.}
        Observet that Vanilla \fedac is indeed less robust to infrequent synchronization and thus worse than the \fedaci. (dataset: \textsc{a9a}, $\lambda=10^{-4}$)}
    \label{fig:instability}
\end{figure}

We include more experiments on various dataset, and more detailed analysis in \cref{additional:expr}.


TITLE: Logarithm Question $ p_1 \cdot p_2 = log(p_1) + log(p_2)? $
QUESTION [0 upvotes]: why are you able to write: $ p_1 \cdot p_2 $ to: $ log(p_1) + log(p_2)  $, where p is the probability of an event occurring.
Also, what would be the benefit of rewriting it in logs?

REPLY [0 votes]: An answer for your second question, you can find here:
http://en.wikipedia.org/wiki/Log_probability


TITLE: A normal to a plane
QUESTION [0 upvotes]: I have been doing one example and there is something unclear. It goes like this:
A plane quadratic net B can be fixed between the points $ (0, 1, 0), (1, 1, 0), (1, 1, −1), (0, 1, −1)$
Now, in the solutions, based on this information it says that: The net B is the subset of the plane y=1 on which 0 ≤ x ≤ 1 and −1 ≤ z ≤ 0. Thus, the constant vector field $\hat N_B$ = (0, 1, 0) is a unit normal field for B.
The thing that I do not understand is: if y=1 is the plane, how is $\hat N_B$ = (0, 1, 0) the normal field to it?
Thanks

REPLY [0 votes]: Rewrite the equation of the plane in vector form as $(x,y,z)\cdot(0,1,0)=0$. This equation says that every point on the plane is orthogonal to $(0,1,0)$, i.e., that $(0,1,0)$ is normal to the plane.


TITLE: Fundamental representation of $\mathfrak{so}(3)$ from the highest weight method
QUESTION [2 upvotes]: Applying the definition of a Lie algebra $\mathfrak{g}$ corresponding to a Lie group $G$ as the tangent space at the identity $T_e(G)$, one can easily construct the Lie algebra $\mathfrak{so}(3)$ corresponding to $SO(3)$ as
$$ \mathfrak{so}(3) = \{ X \in \mathrm{GL}_3(\mathbb{R}):X^T=-X \} $$
i.e. it is the set of $3 \times 3$, real, anti-symmetric matrices. An obvious basis of this space is given by
$$ L_1 = \begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad L_2 = \begin{pmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ -1 & 0 & 0 \end{pmatrix}, \quad L_3= \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & -1 & 0 \end{pmatrix}. $$
Now, if I was to use Cartan's classification and use the highest weight method, I would write down the basis of my fundamental representation as
$$ L_x = \frac{1}{\sqrt{2}} \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \end{pmatrix}, \quad L_y =\frac{1}{\sqrt{2}}\begin{pmatrix}0 & -i & 0 \\ i & 0 & -i \\ 0 & i & 0  \end{pmatrix}, \quad L_z= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0  & -1 \end{pmatrix}. $$
These matrices are neither real nor anti-symmetric. I am aware that the highest weight method produces representations of the Lie algebra, but for the fundamental representation I would expect this to be the trivial case of $ d:\mathfrak{so}(3) \rightarrow \mathrm{GL}(V)$ where $d(X) = X$ for all $X \in \mathfrak{so}(3)$, in other words, I'd expect these matrices to agree, or at least be related by a similarity transformation. As these matrices operate on Euclidean space, there can't be a real similarity transformation to relate these two representations.
My question
Why does the fundamental representation produced by the highest weight method give us a representation that is complex which is neither anti-symmetric or real as the elements of $\mathfrak{so}(3)$ should be?

REPLY [4 votes]: First of all, just for clarity of future readers: the terminology "fundamental representation" is used a bit differently in physics than in mathematics (cf. 1, 2, 3). The representation we are looking at here is the irreducible one of spin $1$ i.e. highest weight $2$; most mathematicians would not count that among what they call "fundamental representations" (for them, the algebra in question has only one fundamental rep, the one with highest weight $1$, which is of dimension $2$). Instead, they might call it the "defining" representation of $\mathfrak{so}_\color{red}3$ (namely, on a $\color{red}3$-dimensional vector space).
The next big obstacle is that further, mathematicians' and physicists' conventions on (in particular compact) Lie algebras differ by multiplying with the imaginary unit $i$. Because I am a mathematician, I'll stick with that notation, differing from the one in the Wikipedia article you cite; meaning that for me,
$$ L_x = \frac{1}{\sqrt{2}} \begin{pmatrix} 0 & i & 0 \\ i & 0 & i \\ 0 & i & 0 \end{pmatrix}, \quad L_y =\frac{1}{\sqrt{2}}\begin{pmatrix}0 & 1 & 0 \\ -1 & 0 & 1 \\ 0 & -1 & 0  \end{pmatrix}, \quad L_z= \begin{pmatrix} i & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0  & -i \end{pmatrix}. $$
Actually, this is necessary, because when you define $\mathfrak{so}_3$ as the anti-symmetric real matrices, you are already following mathematical convention, so we have to stay consistent.
Now even if you grant me that, on first sight your question remains: These are neither real nor anti-symmetric. (They are however anti-hermitian, and make me very happy because comparing some computations I did a long time ago on a bus to this and the higher dimension analogues in that WP article makes me realise that, up to multiplying through with $i$, I did some quantum mechanics without knowing it.)
But now it's conjugation and scaling: With $A:=\dfrac{1}{\sqrt2}\pmatrix{1&0&i\\0&-\sqrt2 i&0\\1&0&-i}$ and $s(M):=AMA^{-1}$, we have
$$s(L_1)= L_x, s(L_3)= L_y, s(L_2)=-L_z . $$
(How to come up with the above base change matrix: Similar to the computations in step 2 in my answer to to Explicit isomorphism between the four dimensional orthogonal Lie algebra and the direct sum of special linear Lie algebras of dimension 3.)
It is not too surprising that the "right" general matrices that come from highest weight theory, i.e. the $L_x, L_y, L_z$ above (in mathematical or physical convention) have complex entries, because the theory relies on complex representations of the complexified Lie algebra, in this case $\mathbb C \otimes \mathfrak{so}_3 \simeq \mathfrak{sl}_2(\mathbb{C})$. What is special for this representation is that there is a similarity transform (the $s$, or rather its inverse, above) which makes it "entirely real" i.e. gives out matrices with all real entries. Of course this would be expected in dimension $3$ because after all we can define $\mathfrak{so}_3$ as a certain set of $3\times3$-matrices with all real entries (and that's what motivated your question). But if instead we had written that Lie algebra as $\mathfrak{su}_2$, that would have been surprising: The standard way one defines $\mathfrak{su}_2$ involves complex $2\times 2$-matrices, and there is no way to transform those into ones with only real entries. If I recall correctly, that parity distinction repeats in higher dimensions, i.e. the irreducible representations in even dimensions (for physicists: half-integer spin) are "truly complex" (or actually, "quaternionic" or "hermitian" or "pseudoreal", depending on what terminology your source likes), whereas the ones in odd dimensions (for physicists: integer spin) can indeed be written with real (or I guess in physicists' notation: purely imaginary) coefficients.


TITLE: How can I evaluate this surface integral
QUESTION [1 upvotes]: Find the value of Surface Integral $ \displaystyle \int \int yzdxdy + xzdydz + xydxdz$ Where the surface is the outer side of the surface formed by the cylinder $x^2 + y^2 =4 $ and the plane $x = 0 , y= 0 ,z = 0$ and  $z=  2$
Now I would like to solve this surface integral with the help of Gauss Divergence theorem, but firstly I want the surface integral in the form $\int F.n ds$
so that I can calculate $\nabla. F$
Can anyone tell me how should I solve this integral ?
Thank you .

REPLY [0 votes]: This is a surface integral of the second kind. A surface integral of the form
$$\iint_R Pdydz + Qdxdz + Rdxdy$$
$$= \iint_S (Pcos\alpha + Qcos\beta + Rcos\gamma)dS$$
Where
$$cos\alpha = \frac{A}{\sqrt{A^2 + B^2 + C^2}}$$
$$cos\beta = \frac{B}{\sqrt{A^2 + B^2 + C^2}}$$
$$cos\gamma = \frac{C}{\sqrt{A^2 + B^2 + C^2}}$$
Where
$$A = \frac{\partial{(y, z)}}{\partial{(u, v)}}$$
$$B = \frac{\partial{(z, x)}}{\partial{(u, v)}}$$
$$C = \frac{\partial{(x, y)}}{\partial{(u, v)}}$$
Based on these formulae we use the following parameter for the cylinder
$$x = 2cos\theta $$
$$y = 2sin\theta $$
$$z = z$$
Using the above formulae we have
$$\sqrt{A^2 + B^2 + C^2} = 2$$
Where
$$A = 2cos\theta$$
$$B = 2sin\theta $$
$$C = 0$$
Substitute x, y, z into the integral we have
$$I = \iint_S (2zsin\theta cos\theta + 2zcos\theta sin \theta)dS$$
Using the change of variables from x, y, z to $\theta$ and z, J factor is 2.
$$= \int_0^{2}\int_0^{\frac{\pi}{2}} 4z sin \theta cos \theta d\theta dz= 4$$


TITLE: Need help in verifying chain condition for applying Zorn's Lemma
QUESTION [0 upvotes]: While self studying abstract algebra from thomas hungerford I have following question in proof of Lemma $6.8$ (Page $222$)

Lemma $6.8.$ Let $A$ be a module over a principal ideal domain $R$ such that $p^{n}A = 0$ and $p^{n-1}A \neq 0$ for some prime $p\in R$ and positive  integer $n$. Let $a$ be an element of $A$ of order $p^n$.
$(i)$ If $A\neq Ra$, then there exists a nonzero $b\in A$ such that $Ra \cap Rb = 0$.
$(ii)$ There is a submodule $C$ of $A$ such that $A = Ra\oplus C$.

Question is in part $2$ of proof $(ii)$

In Case $A\neq Ra$ author defines as $S$ as the set of all submodules $B$ of $A$ such that $Ra\cap B =0$. It is easily proved to be non-empty. Then it is partially ordered by set theoretic inclusion and author says to verify that every chain in $S$ has an upper bound in $S$.

How can I verify that every chain in $S$ has an upper bound in $S$? What upper bound should I take?
I am really confused.
Kindly shed some light on how one should verify upper bound in this sort of questions.

REPLY [4 votes]: Consider $S=\{B\leqslant A\,\mid\,Ra\cap B=0\}$. Let $B_0\subseteq B_1\subseteq B_2\subseteq\dots$ be a chain in $S$. Let $B=\bigcup_{n\in\Bbb N}B_n$. Note that this is the usual first try when dealing with an elementary application of Zorn's Lemma (as already noted by spaceisdarkgreen in the comments).
As any actual calculation in $B$ can already be carried out in some particular $B_{n_0}$ (since the $B_n$ form a chain) for some $n_0\in\Bbb N$, $B_{n_0}$ being a submodule of $A$, we conclude that $B$ is a submodule of $A$ aswell (you might want to make this argument more precise). Now, suppose $Ra\cap B\ne0$ and let $x\in Ra\cap B$ where $x\ne0$. By definition there is some $n_0\in\Bbb N$ such that $x\in B_{n_0}$. But then $Ra\cap B_{n_0}\ne0$ contrary to our assumption. Thus $Ra\cap B=0$ and you can continue by applying Zorn's Lemma as $B$ is an upper bound for the given arbitrary chain.
This sketch is applicable in a lot (not all) similar cases: consider the union of the whole chain and show that this construction yields an upper bound with the same structure.


TITLE: Bound the converging series
QUESTION [0 upvotes]: $ \sum_{n=0}^\infty ((\frac{n^3+2n+1}{n^3+n+1})^{1/3}-1)^{3/4} $
I created a graph and it seems like $ \frac{1}{n} $, but as i know this series should converge. Can't prove it with any usual rules, seems like i have to bound it with another converging function, but couldn't find anything better than $ \frac{1}{n}^{3/4} $, which is obviously don't converge.
$ ((1 +\frac{n}{n^3+n+1})^{1/3}-1)^{3/4} $ is like the most obvious step, but doesn't go that well afterwards. Any throughts?

REPLY [0 votes]: We know that
$$ \lim\limits_{x\to 0} \frac{e^{x/3}-1}{x}  = 1/3 $$
So we can take $x=\ln(\frac{n^3+2n+1}{n^3+n+1}) $. And we know that
$$ \lim\limits_{y\to 0} \frac{\ln(1+y)}{y} = 1 $$
So we can take $y=\frac{n+1}{n^3+n+1}$
Which we can compare further to $\frac{1}{n^2}$. Then by taking the power of $3/4$, we can compare the whole sum to $\frac{1}{n^{3/2}}$, so it's convergent