You might favor location families based on Hamming distance, due to their richness, flexibility, and computational tractability.
Notation and definitions
Recall that in a free finite-dimensional module $V$ with basis $\left(\mathbf{e}_1, \mathbf{e}_2, \ldots, \mathbf{e}_J\right)$, the Hamming distance $\delta_H$ between two vectors $\mathbf{v}=v_1 \mathbf{e}_1 + \cdots + v_J\mathbf{e}_J$ and $\mathbf{w}=w_1 \mathbf{e}_1 + \cdots + w_J\mathbf{e}_J$ is the number of places $i$ where $v_i \ne w_i$.
Given any origin $\mathbf{v}_0\in V$, the Hamming distance partitions $V$ into spheres $S_i(\mathbf{v}_0)$, $i=0, 1, \ldots, J$, where $S_i(\mathbf{v}_0) = \{\mathbf{w}\in V\ |\ \delta_H(\mathbf{w}, \mathbf{v}_0) = i\}$. When the ground ring has $n$ elements, $V$ has $n^J$ elements and $S_i(\mathbf{v})$ has $\binom{J}{i}\left(n-1\right)^i$ elements. (This follows immediately from observing that elements of $S_i(\mathbf{v})$ differ from $\mathbf{v}$ in exactly $i$ places--of which there are $\binom{J}{i}$ possibilities--and that there are, independently, $n-1$ choices of values for each place.)
Affine translation in $V$ acts naturally on its distributions to give location families. Specifically, when $f$ is any distribution on $V$ (which means little more than $f:V\to [0,1]$, $f(\mathbf{v})\ge 0$ for all $\mathbf{v} \in V$, and $\sum_{\mathbf{v}\in V}f(\mathbf{v})=1$) and $\mathbf{w}$ is any element of $V$, then $f^{(\mathbf{w})}$ is also a distribution where
$$f^{(\mathbf{w})}(\mathbf{v}) = f(\mathbf{v}-\mathbf{w})$$
for all $\mathbf{v}\in V$. A location family $\Omega$ of distributions is invariant under this action: $f\in \Omega$ implies $f^{(\mathbf{v})}\in \Omega$ for all $\mathbf{v}\in V$.
Construction
This enables us to define potentially interesting and useful families of distributions by specifying their shapes at one fixed vector $\mathbf{v}$, which for convenience I will take to be $\mathbf{0} = (0,0,\ldots,0)$, and translating these "generating distributions" under the action of $V$ to obtain the full family $\Omega$. To achieve the desired property that $f$ should have comparable values at nearby points, simply require that property of all generating distributions.
To see how this works, let's construct the location family of all distributions that decrease with increasing distance. Because only $J+1$ Hamming distances are possible, consider any decreasing sequence of non-negative real numbers $\mathbf{a}$ = $0 \ne a_0 \ge a_1 \ge \cdots \ge a_J \ge 0$. Set
$$A = \sum_{i=0}^J (n-1)^i\binom{J}{i} a_i$$
and define the function $f_\mathbf{a}:V\to [0,1]$ by
$$f_\mathbf{a}(\mathbf{v}) = \frac{a_{\delta_H(\mathbf{0},\mathbf{v})}}{A}.$$
Then, as is straightforward to check, $f_\mathbf{a}$ is a distribution on $V$. Furthermore, $f_\mathbf{a} = f_{\mathbf{a}'}$ if and only if $\mathbf{a}'$ is a positive multiple of $\mathbf{a}$ (as vectors in $\mathbb{R}^{J+1}$). Thus, if we like, we may standardize $\mathbf{a}$ to $a_0=1$.
Accordingly, this construction gives an explicit parameterization of all such location-invariant distributions that are decreasing with Hamming distance: any such distribution is in the form $f_\mathbf{a}^{(\mathbf{v})}$ for some sequence $\mathbf{a} = 1 \ge a_1 \ge a_2 \ge \cdots \ge a_J \ge 0$ and some vector $\mathbf{v}\in V$.
This parameterization may allow for convenient specification of priors: factor them into a prior on the location $\mathbf{v}$ and a prior on the shape $\mathbf{a}$. (Of course one could consider a larger set of priors where location and shape and not independent, but this would be a more complicated undertaking.)
Generating random values
One way to sample from $f_\mathbf{a}^{(\mathbf{v})}$ is by stages by factoring it into a distribution over the spherical radi and another distribution conditional on each sphere:
Draw an index $i$ from the discrete distribution on $\{0,1,\ldots,J\}$ given by the probabilities $\binom{J}{i}(n-1)^i a_i / A$, where $A$ is defined as before.
The index $i$ corresponds to the set of vectors differing from $\mathbf{v}$ in exactly $i$ places. Therefore, select those $i$ places out of the $\binom{J}{i}$ possible subsets, giving each equal probability. (This is just a sample of $i$ subscripts out of $J$ without replacement.) Let this subset of $i$ places be written $I$.
Draw an element $\mathbf{w}$ by independently selecting a value $w_j$ uniformly from the set of scalars not equal to $v_j$ for all $j\in I$ and otherwise set $w_j=v_j$. Equivalently, create a vector $\mathbf{u}$ by selecting $u_j$ uniformly at random from the nonzero scalars when $j\in I$ and otherwise setting $u_j=0$. Set $\mathbf{w} = \mathbf{v} + \mathbf{u}$.
Step 3 is unnecessary in the binary case.
Example
Here is an R implementation to illustrate.
rHamming <- function(N=1, a=c(1,1,1), n=2, origin) {
# Draw N random values from the distribution f_a^v where the ground ring
# is {0,1,...,n-1} mod n and the vector space has dimension j = length(a)-1.
j <- length(a) - 1
if(missing(origin)) origin <- rep(0, j)
# Draw radii `i` from the marginal distribution of the spherical radii.
f <- sapply(0:j, function(i) (n-1)^i * choose(j,i) * a[i+1])
i <- sample(0:j, N, replace=TRUE, prob=f)
# Helper function: select nonzero elements of 1:(n-1) in exactly i places.
h <- function(i) {
x <- c(sample(1:(n-1), i, replace=TRUE), rep(0, j-i))
sample(x, j, replace=FALSE)
}
# Draw elements from the conditional distribution over the spheres
# and translate them by the origin.
(sapply(i, h) + origin) %% n
}
As an example of its use:
test <- rHamming(10^4, 2^(11:1), origin=rep(1,10))
hist(apply(test, 2, function(x) sum(x != 0)))
This took $0.2$ seconds to draw $10^4$ iid elements from the distribution $f_{\mathbf{a}}^{(\mathbf{v})}$ where $J=10$, $n=2$ (the binary case), $\mathbf{v}=(1,1,\ldots,1)$, and $\mathbf{a}=(2^{11},2^{10},\ldots,2^1)$ is exponentially decreasing.
(This algorithm does not require that $\mathbf{a}$ be decreasing; thus, it will generate random variates from any location family, not just the unimodal ones.)
