Geoscientists often use spatially discretized cellular models of the Earth where data in each grid cell provide independent information about the model parameters of interest at that location. In Bayesian inference this information is given as a set of likelihoods describing the (unnormalized) probability of the parameters, given only the data in each cell. Preexisting information about the model parameters' values and their spatial correlations may be described by a prior probability distribution. The prior, likelihoods, and Bayes' rule together specify a posterior probability distribution that describes the resultant state of information over all model parameters. However, due to the high dimensionality of typical models, the posterior is usually only known up to a multiplicative constant and only at specific, numerically evaluated points in the model space (i.e., it is not known analytically). Markov chain Monte Carlo (MCMC) methods are typically used to produce an ensemble of correlated samples from the posterior. These ensembles are slow to converge in distribution to the posterior; indeed, they may not converge in finite time, and detecting their state of convergence is often impossible in practice. Thus, estimates of the posterior obtained in this way may be biased. We derive a recursive algorithm which samples the posterior exactly, so as to avoid these convergence issues. Its computational cost scales with the size of the parameters' sample space, the prior's spatial range of dependency, and the shortest edge dimension of the grid. We develop an approximation to the algorithm such that it may be used on large 2‐D (and potentially 3‐D) model grids. We apply it to synthetic seismic attribute data and obtain results which compare favorably to the results of MCMC (Gibbs) sampling—which exhibits convergence problems.