How can I create a map of variance for IDW predictions? Unlike with kriging, predicting spatial variation using the IDW function from gstat returns only predicted values, not estimates of variation.  
library(sp)
data("meuse")
coordinates(meuse) = ~x+y
data("meuse.grid")
gridded(meuse.grid) = ~x+y
m.idw <- gstat::idw(zinc~1, meuse, newdata=meuse.grid, idp=2)

There is a column for variation, but it is empty.  Is this a fundamental mathematical limitation of the IDW approach, or is it just not implemented in gstat, or is there something else wrong in my approach?
 A: I will supply two answers: first, this is a limitation of IDW; but second, there are ways to overcome it by making assumptions.

IDW (inverse distance weighting) is a procedure--best used in an exploratory mode--for predicting unobserved values of a function defined on a metric space such as the Euclidean plane or sphere.  Specifically, suppose you have observed $y_i = f(x_i)$ at a set of distinct points $x_i\in \mathcal{X}$ where $\mathcal X$ has a metric $\delta$ and you wish to predict the values of $f$ at other points $x^\prime_j\in \mathcal X$ (which often are chosen to support the creation of graphs, contour maps, images, etc. of $f$).  
Carrying out IDW requires two choices of the analyst (both of which are often made arbitrarily by accepting software defaults): (1) selecting a power $p\ge 1$ and (2) defining a neighborhood search procedure.  The latter is a way of selecting observation points near any given point $x.$  For our purposes the details don't matter, so let's just write $\mathcal{N}(x) \subset \{x_1,x_2,\ldots, x_n\}$ for the neighborhood of $x\in\mathcal{M}.$
Given these choices, IDW is just a recipe: from the distances $\delta(x,x_i)$ between $x$ and its neighboring observations $x_i\in\mathcal{N}(x),$ construct relative weights 
$$w_i(x) = \delta(x,x_i)^{-p} \gt 0,$$
normalize them to sum to unity by setting
$$\omega_i(x) = \frac{w_i(x)}{\sum_{x_i\in\mathcal{N}(x)} w_i(x)},$$
and predict the function's value as a weighted combination of the neighboring observations,
$$\hat y(x) =\sum_{x_i\in\mathcal{N}(x)} \omega_i(x)\, y_i.$$

Because there is no probability model involved, there is no general way to estimate a variance or standard error of estimate.


One could use various cross-validation tactics to estimate prediction errors.  These all implicitly assume $f$ has various restrictive properties concerning how rapidly it can change from point to point.  Rather than dodging considerations of what is actually being assumed about $f,$ I have found it useful at times to do a quick-and-dirty geostatistical analysis of the data, which (although it also makes assumptions) can reveal something about the properties of $f.$
This means supposing $f$ is a realization of some kind of stochastic process $F$ ("random field") on $\mathcal{M},$ so that $(y_1,y_2,\ldots,y_n)$ is a single observation of the $n$-variate random variable $(F(x_1), F(x_2), \ldots, F(x_n));$ and estimating its covariance function
$$C(x, x^\prime) = \operatorname{Cov}(F(x), F(x^\prime)).$$
This is usually accomplished by adopting some assumptions about stationarity and ergodicity so that the data can be used to estimate $C$ (in the usual way, likely through some form of variogram estimation).  At other times, you might use estimates of $C$ based on related data.

With such a covariance function in hand it's a relatively simple matter to compute the variance of any predicted value as $$\operatorname{Var}(\hat{y}(x)) = \operatorname{Var}\left(\sum_{x_i\in\mathcal{N}(x)} \omega_i(x)\, y_i\right) = \sum_{i,j\in\mathcal{N}(x)} \omega_i(x)\omega_j(x)\,C(x_i,x_j).$$

Its square root serves as an estimated standard error of prediction.
Note that what makes the foregoing formula possible is that the weights $\omega_i(x)$ depend only on the locations of the observations, which are considered determinate (non-random), and not on the observations themselves.

Among other applications, this approach affords a principled way to adjust the otherwise arbitrary choices of power $p$ and neighborhood search $\mathcal{N}$ made at the outset: if one set of choices consistently results in smaller standard errors (in regions of $\mathcal M$ of particular interest) than another set, then the former is to be preferred.  Note that making such a comparison does not require excellent estimation of $C$; in particular, the hardest part of estimating $C$ usually is finding an appropriate scale factor--but that's irrelevant if one is only comparing two predictions based on the same $C.$  Usually, determining the right range of the covariance function is what matters most.  That's why this approach usually is not onerous and can produce useful information.
A: As others already pointed out, IDW is not based on a theory that could provide variance maps. On the other hand, it is easy to arbitrarily choose a metric, and use it to quantify how far from the observations the estimates are. This is what I have implemented in GeoStats.jl:
using GeoStats
using InverseDistanceWeighting
using Plots

sdata   = PointSetData(Dict(:z => [1.,0.,1.]), [25. 50. 75.; 25. 75. 50.])
sdomain = RegularGrid{Floa64}(100, 100)

problem = EstimationProblem(sdata, sdomain, :z)

solver = InvDistWeight()

solution = solve(problem, solver)

contourf(solution)


And this is how it looks like if you choose a Chebyshev metric:
solver = InvDistWeight(:z => (distance=Chebyshev(),))

solution = solve(problem, solver)

contourf(solution)


