One option may be to split the original data into two subsets: one that will be used in interpolating values and one that will be used to validate the interpolation results. The error is then estimated by comparing interpolated values at the validation point locations with the actual validation point values. Note that the appropriateness of this approach is largely driven by the sample point density and distribution vis-à-vis the type and scale of the underlying process you are attempting to model.
Edit: This is an expansion of the original answer following @whuber’s comments.
As noted by @whuber, one disadvantage with the aforementioned technique is the degradation in the interpolation’s quality when a subset of the sampled sites are removed. A solution to this problem is described in Maciej Timczak’s 1998 paper. The author applies a cross-validation technique to estimate the optimal interpolation parameters. He then uses a jackknife method to estimate the predicted value and uncertainty at an unsampled site $Z_J$. A brief summary of the techniques described in the paper with a simple example in R follows.
In the cross-validation (leave-one-out) method, one data point $s_i$ is removed from the point data set $S$ and its interpolated value is computed using all other $(n-1)$ points of $S$. The interpolated value is then compared with the actual value $s_i$. This process is repeated with all other data points from $S$. The performance of the interpolator is evaluated via the root-mean of squared residuals (RMSE). The RMSE can be computed for different interpolation parameters (or different interpolators) then compared. The interpolator with the lowest RMSE is usually desired.
$$RMSE=\sqrt{\frac{\sum_{i=1}^n (Z_{i(int)} -Z_{i})^2}{n}}$$
Once the interpolation technique is chosen, the jackknife technique is used to estimate the unknown value $Z_j$ at an unsampled location $s_j$ along with its confidence interval. The method, as implemented by the author, involves first interpolating the unknown value $Z_j$ at unsampled location $s_j$ using all sample points from $S$, then interpolating a pseudo-value $Z_i^*$ using $(n-1)$ points from the sample dataset $S$ (i.e. $Z_i^*$ is computed once for each omitted point sample $s_i$):
$$Z_{i}^* = n Z_{all} - (n-1) Z_{-i}$$
where i = 1,2,..,n
The jackknifed estimator of $Z_j$ at $s_j$ is computed by averaging all of the pseudo-values following:
$$Z_{J} =\frac{ \sum_{i=1}^{n} Z_{i}^*}{n}$$
A confidence interval is then computed as follows:
$$\sigma_J=\sqrt{\frac{1}{n(n-1)}\sum_{i=1}^{n}(Z_i^*-Z_J)^2}$$
The estimated value is therefore $Z_j \pm- t_{(\alpha/2,n-1)}\sigma_J$
A simple R example follows:
X <- c(100.0,19.4,9.0,64.4,39.4,50.7,99.0,44.4,82.5,55.9,
56.2,54.0,14.9,54.8,35.5,34.6,15.2,32.0,23.8,87.4,
49.4,77.9,63.9,14.8,5.9,45.3,95.6,10.3,59.5,47.2,
26.7,46.5,41.3,62.9,34.2,3.7,57.7,78.5,73.1,28.3,
13.1,49.4,24.2,99.2,76.3,93.2,71.6,28.8,49.4,94.0,
84.4,0.0,90.3,48.4,44.8,5.1,29.8,27.7,93.8,25.6)
Y <- c(0.0,1.9,3.2,6.0,12.4,13.3,13.7,15.3,15.6,
16.0,18.0,22.3,22.9,23.0,24.3,26.3,26.6,27.4,
31.6,33.1,33.8,35.0,35.2,42.0,44.9,45.3,45.8,
48.8,50.5,58.3,60.2,60.8,60.8,61.5,64.4,64.6,
65.5,69.2,69.3,69.4,71.2,73.3,78.5,80.1,83.6,
84.7,84.7,91.1,92.1,92.7,93.0,93.7,93.8,95.1,
96.0,96.4,96.7,97.6,99.8,100.0)
Z <- c(209,478,424,817,866,720,327,833,731,1488,562,868,318,496,488,
1146,369,735,593,778,771,304,538,669,368,474,391,346,872,556,
348,765,779,809,357,720,416,544,338,560,455,555,340,307,589,
280,745,452,1116,442,659,343,385,655,828,490,425,665,276,333)
library(akima)
n <- 60 # Number of points in S
#
# Cross validation (leave-one-out) method to compare interpolators
#
P.spl <- vector()
sum.dif2 <- numeric()
for (i in 1:n){
P.spl[i] <- interpp(X[-i],Y[-i],Z[-i],X[i],Y[i],linear=F,extrap=T)$z
sum.dif2 <- (P.spl[i] - Z[i])^2
}
rmse = sqrt(sum.dif2 / n)
rmse
#
# Jackknife to estimate the uncertainty in the interpolated value
P.spl <- vector()
Zi <- numeric()
jx <- 40 # X coordinate of parameter Zj to be estimated at sj
jy <- 40 # Y coordinate of parameter Zj to be estimated at sj
Zall <- interpp(X,Y,Z,jx,jy,linear=T)$z # Interp value from all si
for (i in 1:n){
Z1 <- interpp(X[-i],Y[-i],Z[-i],jx,jy,linear=T)$z # Interp value from si-1
# Calculated pseudo-value Z at j
Zi[i] <- n * Zall - (n-1) * Z1
}
#
# Jackknifed estimator of parameter Z at location j
#
Zj <- sum(Zi) / n # Estimated value Zj
# Estimated standard error
sig.j <- sqrt(1/(n*(n-1)) * sum((Zi-Zj)^2))
# The confidence interval on Zj
alpha <- 1 - .05 / 2
t.value <- qt(alpha,n-1)
ci <- t.value * sig.j # Confidence interval for Zj