{'fraction_non_alphanumeric': 0.04052737287967069, 'fraction_numerical': 0.011056535486308588, 'mean_word_length': 4.15615707987286, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Prior specification for binary Markov mesh models

We propose prior distributions for all parts of the specification of a Markov mesh model. In the formulation we define priors for the sequential neighborhood, for the parametric form of the conditional distributions and for the parameter values. By simulating from the resulting posterior distribution when conditioning on an observed scene, we thereby obtain an automatic model selection procedure for Markov mesh models. To sample from such a posterior distribution, we construct a reversible jump Markov chain Monte Carlo algorithm (RJMCMC). We demonstrate the usefulness of our prior formulation and the limitations of our RJMCMC algorithm in two examples.


Introduction
Discrete Markov random fields (MRFs) and Markov mesh models (MMMs) defined on rectangular lattices are popular model classes in spatial statistics, see for example Kindermann and Snell (1980) and Hurn et al. (2003) for MRFs, and Abend et al. (1965) and Cressie and Davidson (1998) for MMMs. Discrete MRFs are frequently used to model available prior information about an unobserved scene x of a discrete variable. This prior is then combined with a likelihood function describing the relation between x and some observed data y into a posterior distribution, and this posterior is the basis for making inference about x.
When specifying the MRF prior, the most frequent approach is to fix the neighborhood and parametric model structures and also to specify the values of the model parameters a priori. However, some authors have also explored a more fully Bayesian approach (Heikkinen and Högmander, 1994;Higdon et al., 1997;Friel et al., 2009;Everitt, 2012). In particular, Arnesen and Tjelmeland (2017) formulate a prior for all parts of the MRF prior and demonstrate how MCMC sampling from the corresponding posterior distribution when conditioning on an observed scene produces MRFs that give realizations with a similar spatial structure as present in the scene used to define the posterior.
The class of Markov mesh models, and the partially ordered Markov model (POMM) generalization (Cressie and Davidson, 1998) of this class, is much less used in the literature. We think the main reason for this is that it is much harder to manually choose an MMM than an MRF that reflects given prior information.
It is neither an easy task to specify an MRF that is consistent with given prior information, but except for boundary effects it is for an MRF easy to ensure that the field is stationary. This is an important practical argument when completely specifying the prior a priori, but it is not so important when a fully Bayesian model is adopted as in Arnesen and Tjelmeland (2017). It should also be noted that MRFs contain a computationally intractable normalizing constant which severely limits the practicability of MRFs in a fully Bayesian context, see for example the discussion in Friel et al. (2009). In contrast, the normalizing constant for an MMM is explicitly given in an easy to compute form. Also for this reason an MMM is much better suited as a prior than an MRF when adopting the fully Bayesian approach.
Our goal in the present article is to formulate a fully Bayesian MMM. In particular, we would like the hyper-prior to include distributions for the neighborhood structure, for the interaction structure of the conditional distributions defining the MMM, and for the parameter values. We should thereby obtain a flexible prior that is able to adapt to a wide variety of scenes. To specify the MMM hyperprior, we adapt the general strategy used in Arnesen and Tjelmeland (2017) for the MRF to our MMM situation. Given such a Bayesian model, we also want to formulate an MCMC algorithm to simulate from the resulting posterior distribution conditioned on an observed scene. It should thereby be possible to learn both the form of the parametric model and the values of the model parameters from an observed scene. For simplicity we here limit our attention to binary MMMs, but our approach can quite easily be generalized to a situation where each node has more than two possible values.
The remainder of this article is organized as follows. In Section 2 we introduce most of the notations we use for defining our Bayesian Markov mesh model, and in particular discuss pseudo-Boolean functions. In Section 3 we use this to formulate the Markov mesh model class. In Section 4 we construct our prior distribution, and in Section 5 we formulate proposal distributions that we use in a reversible jump Markov chain Monte Carlo algorithm to simulate from the corresponding posterior when conditioning on an observed scene. In Section 6 we present two simulation examples and lastly we give some closing remarks in Section 7.

Preliminaries
In this section we first introduce the notation we use to represent a rectangular lattice, the variables associated to this lattice and some quantities we use to formulate our Markov mesh model defined on this lattice. Thereafter, we define the class of pseudo-Boolean functions and explain how a pseudo-Boolean function can be used to represent a conditional distribution for binary variables.

Notation
Consider a rectangular m × n lattice. Let v = (i, j) denote a node in this lattice, where i and j specify the vertical and horizontal positions of the node in the lattice, respectively. We let i = 1 be at the top of the lattice and i = m at the bottom, and j = 1 and j = n are at the left and right ends of the lattice, respectively. We use lowercase Greek letters to denote sets of nodes, and in particular we let χ = {(i, j) : i = 1, . . . , m, j = 1, . . . , n} be the set of all nodes in the lattice. Occasionally we also consider an infinite lattice Z 2 , where Z = {0, ±1, ±2, . . .} is the set of all integers, and we use v = (i, j) ∈ Z 2 also to denote a node in such an infinite lattice. We use λ, λ ⋆ ⊆ Z 2 to denote arbitrary sets of nodes. To translate a node (i, j) ∈ Z 2 by an amount (k, l) ∈ χ, we adopt the notation (i, j) ⊕ (k, l) = (i + k, j + l). (1) One should note that even if (i, j) ∈ χ, (i, j) ⊕ (k, l) may fall outside the finite m × n lattice. To translate all nodes in a set λ ⊆ Z 2 by the same amount (k, l) ∈ χ, we write λ ⊕ (k, l) = {(i, j) ⊕ (k, l) : (i, j) ∈ λ}. (2) To denote sets of subsets of nodes, we use uppercase Greek letters, and in particular, we let Ω(χ) = {λ : λ ⊆ χ} denote the set of all subsets of χ, often called the power set of χ. One should note that Ω(χ) in particular includes the empty set and χ itself. We use Λ, Λ ⋆ ⊆ Ω(χ) to denote arbitrary sets of subsets of nodes.
To define a Markov mesh model one must, for each node v = (i, j), define a socalled predecessor set and a sequential neighborhood. After numbering the nodes in the lattice from one to mn in the lexicographical order, we let the predecessor set of a node (i, j) consist of all nodes with a lower number than the number of see the illustration in Figure 1(a). We let ν v = ν (i,j) ⊆ ρ (i,j) denote the sequential neighborhood for node v = (i, j) ∈ χ as illustrated in Figure 1(b). In Section 3 we consider a Markov mesh model where all the sequential neighborhoods are defined by a translation of a template sequential neighborhood τ . The τ can be thought of as the sequential neighborhood of node (0, 0) in the infinite lattice. More precisely, τ is required to include a finite number of elements and where Z − = {−1, −2, . . .} is the set of all negative integers. The sequential neighborhood of a node v ∈ χ is then defined as As illustrated in Figure 2, sequential neighborhoods for all nodes sufficiently far away from the lattice borders then have the same form, whereas nodes close to the borders have fewer sequential neighbors. One can note that with this construction one always has ν (1,1) = ∅.
To each node v = (i, j) ∈ χ, we also associate a corresponding binary variable which we denote by The collection of all these binary variables we denote by x = (x v ; x ∈ χ) and we let x λ = (x v ; v ∈ λ) represent the collection of variables associated to the nodes in a set λ ⊆ χ. In particular, x νv is the collection of variables associated to the sequential neighborhood of node v. If x v = 1 we say node v is on, and if x v = 0 we say the node is off. We let ξ(x) ⊆ χ denote the set of all nodes that are on, i.e.
In particular, the set of nodes in the sequential neighborhood of node v that is on is then, using (5) and that ξ(x) ⊆ χ, In the next section, we define the class of pseudo-Boolean functions which we in Section 3 use to define the class of Markov mesh models.

Pseudo-Boolean functions
When defining pseudo-Boolean functions, we reuse some of the symbols introduced when discussing concepts related to the rectangular m × n lattice above.
In particular, we define a pseudo-Boolean function with respect to some finite set, denoted by τ . In the definition, this τ has no relation to the template sequential neighborhood τ introduced above. However, when applying a pseudo-Boolean function to represent the conditional distribution of x v for a node v ∈ χ given the values of the nodes in the sequential neighborhood ν v , the set τ used to define a pseudo-Boolean function is equal to the template sequential neighborhood τ .
In particular, the elements of τ is then the nodes in the lattice χ, and therefore we use λ and λ ⋆ to represent subsets of τ also when discussing pseudo-Boolean functions in general.
A pseudo-Boolean function θ(·) defined on a finite set τ is a function that associates a real value to each subset of τ , i.e.
where Ω(τ ) is the power set of τ . Thereby, for any λ ⊆ τ the value of the pseudo-Boolean function is θ(λ). Equivalently, one may think of a pseudo-Boolean function as a function that associates a real value to each vector z ∈ {0, 1} |τ | , where |τ | is the number of elements in the set τ . To see the correspondence, one should set an element in z equal to one if and only if the corresponding element in τ is in the set λ. This last formulation of pseudo-Boolean functions is the more popular one, see for example Hammer and Holzman (1992) and Grabisch et al. (2000), but in the present article we adopt the formulation in (8) as this gives simpler expressions when formulating our Markov mesh model in Section 3 and the corresponding prior distribution in Section 4. Hammer and Rudeanu (1968) show that any pseudo-Boolean function can be uniquely represented by a collection of interaction parameters (β(λ), λ ∈ Ω(τ )) by the relation The corresponding inverse relation is given by The one-to-one relation in (9) and (10) is known as Moebious inversion, see for example Lauritzen (1996).
If one or more of the interaction parameters β(λ) are restricted to be zero, a reduced representation of the pseudo-Boolean function can be defined. For some Λ ⊆ Ω(τ ) assume now that one restricts β(λ) = 0 for all λ ∈ Λ. One can then represent the pseudo-Boolean function θ(·) by the interaction parameters {β(λ), λ ∈ Λ}, and the relation in (9) becomes where Ω(λ) is the power set of λ. We then say that θ(·) is represented on Λ.
Moreover, we say that the set Λ is dense if for all λ ∈ Λ, all subsets of λ is also included in Λ, and that the template sequential neighborhood τ is minimal for Λ if all nodes v ∈ τ are included in at least one of the elements of Λ. One should note that if Λ is dense and τ is minimal for Λ then there is a one-to-one relation between the elements in τ and the sets λ ∈ Λ which contains only one node, Throughout this paper, we restrict attention to pseudo-Boolean functions that are represented on a Λ that is dense and the template sequential neighborhood τ that is minimal for this Λ. A λ ∈ Ω(τ ) we term an interaction, we say the interaction is active if λ ∈ Λ and otherwise we say it is inactive. The Λ is thereby the set of active interactions.
As also discussed in Austad and Tjelmeland (2017), the set of active interactions Λ can be visualized by a directed acyclic graph (DAG), where we have one to represent the node (0, 0), whereas each node (i, j) ∈ λ for each λ ∈ Λ is represented by a placed in position (i, j) relative to ⊠.
vertex for each active interaction λ ∈ Λ and a vertex λ ∈ Λ is a child of another vertex λ ⋆ ∈ Λ if and only if λ = λ ⋆ ∪{v} for some v ∈ τ \λ ⋆ . Figure  This τ can be used to define the sequential neighborhoods for nodes in a rectangular lattice as discussed in Section 2.1. In the vertices of the DAG shown in the figure, node (0, 0) is represented by the symbol ⊠, whereas each of the nodes in λ ∈ Λ is represented by the symbol . Thinking of τ as a finite set of nodes in a lattice, the position of the representing node (i, j) ∈ λ is placed at position (i, j) relative to ⊠.

Markov mesh model
In this section we formulate a homogeneous binary Markov mesh model (Abend et al., 1965) for a rectangular m × n lattice. We adopt the notation introduced in Section 2, so in particular χ denotes the set of all nodes in the m × n lattice and is the collection of the binary variables associated to χ. In a Markov mesh model the distribution of x is expressed as where f (x v |x ρv ) is the conditional distribution for x v given the values of the variables in the predecessor nodes. Moreover, one assumes the Markov property i.e. the conditional distribution of x v given the values in all predecessors of v only depends on the values in the nodes in the sequential neighborhood of v. As discussed in Section 2.1, we assume the sequential neighborhoods ν v , v ∈ χ to be defined as translations of a template sequential neighborhood τ as specified in (4) and (5). Using the result in (7), the conditional distribution f (x v |x νv ) can then be uniquely represented by a pseudo-Boolean function θ v (λ), λ ⊆ τ by the In general, one may have one pseudo-Boolean function θ v (λ) for each v ∈ χ, but in the following we limit the attention to homogeneous models, so we require all θ v (·), v ∈ χ to be equal. We let θ(·) denote this common pseudo-Boolean function, i.e. θ v (λ) = θ(λ) for all λ ⊆ τ and v ∈ χ and, without loss of generality, we assume θ(·) to have a dense representation on a set Λ ⊆ Ω(τ ) and τ to be minimal for Λ. Thus, the distribution of our homogeneous binary Markov mesh model is Assuming, as we do, the Markov mesh model to be homogeneous is convenient in that we do not need to specify a separate pseudo-Boolean function for each node v ∈ χ, and it is also statistically favorable as it limits the number of parameters in the model. However, one should note that this choice implies that for a node v ∈ χ close to the boundary of the lattice so that the set (τ ⊕ v) \ χ is non-empty, the conditional distribution f (x v |x νv ) is as if the nodes (for the infinite lattice) in the translation of τ that fall outside the lattice χ are all zero. Thus, even if the model is homogeneous it is not stationary, and in particular one should expect strong edge effects since we in some sense are conditioning on everything outside the lattice χ to be zero. When estimating or fitting the model to an observed scene, it is crucial to take this edge effect into account.
Having defined our class of homogeneous Markov mesh models as above, a model is specified by the template sequential neighborhood τ , the set of active interactions Λ ⊆ Ω(τ ) on which the pseudo-Boolean function θ(·) is represented, and the parameter values {θ(λ) : λ ∈ Λ}. Thus, to adopt a fully Bayesian approach, we need to formulate prior distributions for τ , Λ and {θ(λ) : λ ∈ Λ}, and this is the focus of the next section.

Prior distribution
When constructing our prior distribution for the template sequential neighborhood τ , the set of active interactions Λ and the parameter values {θ(λ) : λ ∈ Λ}, we have two properties in mind. Firstly, the prior should be vague so that the Markov mesh model manages to adapt to a large variety of scenes. To obtain this, the number of elements in τ should be allowed to be reasonably large and higherorder interactions should be allowed in the model. Secondly, to avoid overfitting, the prior should favor parsimonious Markov mesh models, and in particular this implies that the highest prior probabilities should be assigned to models with just a few higher-order interactions.
We define the prior as a product of three factors where f (τ ) is a prior for the template sequential neighborhood τ , f (Λ|τ ) is a prior for the set of active interactions Λ when τ is given, and is a prior for the parameter values given τ and Λ. In the following we discuss each of these factors in turn.

Prior for the template sequential neighborhood τ
We restrict the template sequential neighborhood to be a subset of a given finite set τ 0 ⊂ ψ, where ψ is defined in (4). The τ 0 can be though of as a set of possible sequential neighbors for node (0, 0). To get a flexible prior it is important that the number of elements in τ 0 is not too small, and it is natural to let τ 0 include nodes close to (0, 0). For example, one may let ψ include all nodes that are inside the circle centered at (0, 0) with some specified radius r.
In the examples discussed in Section 6 we use this with r = 5, see the illustration in Figure 4.
Given the set τ 0 we specify the prior for τ ⊆ τ 0 by first choosing a prior distribution for the number of elements in τ , and thereafter a prior for τ given the number of elements in τ . Letting n τ = |τ | denote the number of elements in τ we thereby have For simplicity we choose both f (n τ ) and f (τ |n τ ) to be uniform distributions. The possible values for n τ are all integers from 0 to |τ 0 |, so we get Moreover, having chosen τ to be uniform given n τ = |τ |, we get where the binomial coefficient in the numerator is the number of possible sets τ 's with n τ elements.
One should note that our choice of the two uniforms above is very different from adopting a uniform prior for τ directly. A uniform prior on τ would have resulted in very high a priori probabilities for n τ being close to |τ 0 |/2 and very small a priori probabilities for values of n τ close to zero, which is clearly not desirable.
One can easily construct other reasonable priors for τ than the one defined above. For example, one could want to build into the prior f (τ |n τ ) that nodes close to (0, 0) are more likely to be in τ than nodes further away. Recalling that we want to simulate from a corresponding posterior distribution by a reversible jump Markov chain Monte Carlo algorithm (RJMCMC) (Green, 1995), the challenge is to formulate a prior with this property so that we are able to compute the (normalized) probability f (τ |n τ ), as this is needed to evaluate the Metropolis-Hastings acceptance probability. For the data sets discussed in Section 6, we have also tried a prior f (τ |n τ ) in which we split the nodes in τ 0 into two or three zones dependent on their distances from (0, 0) and have a different prior probability for a node to be in τ dependent on which zone it is in. As long as the number of zones is reasonably small, it is then possible to compute the normalizing constant of f (τ |n τ ) efficiently. However, in our simulation examples this gave essentially the same posterior results as the very simple double uniform prior specified above.

Prior for the set of active interactions Λ
To specify a prior for the set of active interactions Λ, we first split Λ into several subsets dependent on how many nodes an element λ ∈ Λ contains. More precisely, Thus, Ω k (τ ) contains all k'th order interactions, and Λ k ⊆ Ω k (τ ) is the set of all k'th order active interactions. As we have assumed τ to be minimal for Λ, τ is uniquely specifying Λ 1 = {λ ∈ Λ : |λ| = 1}, see the discussion in Section 2.2 and in particular (12). Moreover, we restrict ∅ always to be active, i.e. ∅ ∈ Λ with probability one, which implies that we force the pseudo-Boolean function θ(·) always to include a constant term. As we have already assumed Λ to be dense and τ to be minimal for Λ this is only an extra restriction when τ = ∅. Thus, for given τ the sets Λ 0 and Λ 1 are known, so to formulate a prior for Λ we only need to define a prior for Λ k , k = 2, . . . , |τ |. We assume a Markov property for these sets in that Thus, to choose a prior f (Λ|τ ) we only need to formulate f (Λ k |Λ k−1 ), and to do so we adopt the same strategy for all values of k. In the specification process of f (Λ k |λ k−1 ) we should remember that we have already restricted Λ to be dense, so the chosen prior needs to be consistent with this. For a given Λ k−1 , an interaction are active. We let Π k denote this set of possible active k'th order interactions, i.e. we must have We assume each interaction λ ∈ Π k to be active with some probability p k , independently of each other, and get One should note that if Λ k−1 = ∅ one gets Π k = ∅ and thereby also f (Λ k = ∅|Λ k−1 ) = 1.
The probabilities p k , k = 2, . . . , |τ | should be chosen to get a reasonable number of higher-order active interactions. To obtain a parsimonious model, one need to adopt a small value for p k if the number of elements in Π k is large, but to favor a model to include some higher-order interactions, the value of p k can be large when the number of elements in Π k is small. We choose where p ⋆ ∈ (0, 1) is a hyper-parameter to be specified. One should note that this choice in particular ensures the expected number of active k'th order interactions to be smaller than |Λ k−1 |.

Prior for the parameter values {θ(λ) : λ ∈ Λ}
Given τ and the set of active interactions Λ, the set of model parameters for which we need to formulate a prior is {θ(λ) : λ ∈ Λ}. From the model assumptions in (15) and (16), we have that each θ(λ), λ ∈ Λ have a one-to-one correspondence with the conditional probability Since the θ(λ)'s define probabilities conditioning on different values for x ρv , we find it reasonable, unless particular prior information is available and suggests otherwise, to assume the θ(λ), λ ∈ Λ to be independent. In the following we adopt this independence assumption. Moreover, as we do not have a particular class of scenes in mind but want the prior to be reasonable for a wide variety of scenes, we adopt the same prior density for all parameters θ(λ), λ ∈ Λ.
To formulate a reasonable and vague prior for θ(λ), we use the one-to-one correspondence between θ(λ) and the probability p(λ). The interpretation for p(λ) is much simpler than that of θ(λ), so our strategy is first to choose a prior for p(λ) and from this derive the corresponding prior for θ(λ). As we do not have a particular class of scenes in mind but want our prior to be reasonable for a wide variety of scenes, we find it most natural to adopt a uniform prior on [0, 1] for p(λ). However, as previously mentioned we want to explore a corresponding posterior distribution by running a reversible jump Metropolis-Hastings algorithm, and in particular we want to use adaptive rejection sampling (Gilks, 1992) to update θ(λ). For this to work, the full conditional for θ(λ) needs to be log-concave.
Adopting the uniform on [0, 1] prior for p(λ) the resulting posterior full conditional becomes log-concave, but the second derivative of the log full conditional converges to zero when θ(λ) goes to plus or minus infinity. As this may generate numerical problems when running the adaptive rejection sampling algorithm, we adopt a prior for p(λ) slightly modified relative to the uniform and obtain a posterior distribution where the second derivative of the log full conditional for θ(λ) converges to a value strictly less than zero. More precisely, we adopt the following prior for θ(λ), where the first factor is the prior resulting from assuming p(λ) to be uniform, the second factor is the modification we adopt to avoid numerical problems when running the adaptive rejection sampling algorithm, and σ > 0 is a hyper-parameter to be specified. The resulting priors for p(λ) and θ(λ) when σ = 10 are shown in Figure 5. We see that the prior for p(λ) is close to the uniform. One can also note that f (θ(λ)|Λ) have heavier tails than a normal distribution with the same variance. One should note that the normalizing constant in (27)  given in (27).
Having specified priors for τ , Λ and {θ(λ) : λ ∈ Λ} we formulate in the next section a reversible jump Metropolis-Hastings algorithm for simulating from the corresponding posterior when a scene x is observed.

Simulation algorithm
In the section we assume we have observed a complete scene x = (x v ; v ∈ χ) and assume this to be a realization from the Markov mesh model defined in Section 3. We adopt the prior defined in Section 4 and want to explore the resulting posterior distribution by a reversible jump Markov chain Monte Carlo algorithm (RJMCMC), see Green (1995). We combine two types of updates. In the first update class, we keep τ and Λ unchanged and update the parameter vector {θ(λ) : λ ∈ Λ} by a Gibbs step along a direction sampled uniformly at random. In the second update class, we propose a trans-dimensional move by adding an inactive interaction to Λ or removing an active interaction from Λ, and proposing corresponding changes for the parameter vector {θ(λ) : λ ∈ Λ}.
It is clearly of interest to consider also the resulting posterior distribution when parts of the scene x is unobserved or when x is an unobserved latent field.
The former is of interest if one wants to reduce the boundary effects of the Markov mesh model by letting x include an unobserved boundary around the observed area, and the latter is a common situation in image analysis applications. However, to simplify the discussion of the simulation algorithm in this section, we assume the complete scene x to be observed. In Section 6, where we present a number of simulation examples, we describe how to adapt the simulation algorithm to situation in which a part of x is unobserved.
In the following we describe each of the two update types in turn, starting with the Gibbs update for the parameter values. We only discuss the proposal distribution, as the acceptance probabilities is then given by standard formulas. by sampling ∆(λ) from a standard normal distribution, independently for each λ ∈ Λ. We then set where α ∈ R is sampled from the full conditional As α is sampled from its full conditional, this is a Gibbs update and the Metropolis-Hastings acceptance probability is one. The full conditional (31) for α is not of a standard form, but in Appendix A we show that it is log-concave, so to generate samples from it we adopt the adaptive rejection sampling algorithm of Gilks (1992).

Updating the set of active interactions
Again let again τ , Λ and {θ(λ) : λ ∈ Λ} be the current state. In this update we modify Λ, and possibly also τ , by adding an inactive interaction to Λ or by removing an active interaction from Λ. We let τ ⋆ and Λ ⋆ denote the potential new values for τ and Λ, respectively. With a change in Λ, the number of parameter values {θ(λ) : λ ∈ Λ} is also changed, and to try to obtain a high acceptance rate, we in fact propose a change also in some of the parameter values that are in both the current and potential new states. We let {θ ⋆ (λ) : λ ∈ Λ ⋆ } denote the set of potential parameter values.
To generate τ ⋆ , Λ ⋆ and {θ ⋆ (λ) : λ ∈ Λ ⋆ }, we first draw at random whether to add an inactive interaction to Λ or to remove an active interaction from Λ. In the following we specify in turn our procedures for proposing to remove and add an interaction.
By setting the algorithmic tuning parameter ν = 0, we draw Λ ⋆ uniformly at random from the elements in Λ r . With a larger value for ν, we get higher probability for proposing to remove an interaction λ ⋆ that gives a small change in the pseudo-Boolean function. If it should happen that Λ r = ∅, we simply propose an unchanged state. Assuming we have sampled a λ ⋆ to remove, we have two possibilities. If λ ⋆ is a higher-order interaction the sequential neighborhood is unchanged, i.e. τ ⋆ = τ , whereas if λ ⋆ is a first-order interaction the sequential neighborhood is reduced to τ ⋆ = τ \ λ ⋆ .
Having decided τ ⋆ and Λ ⋆ , the next step is to specify the potential new parameter values {θ ⋆ (λ) : λ ∈ Λ ⋆ }. To understand our procedure for doing this, one should remember that there is a one-to-one relation between the current parameter values {θ(λ) : λ ∈ Λ} and a set of current interaction parameters where the relation is given by (9) and (10) with respect to {θ ⋆ (λ) : λ ∈ Ω(τ 0 )}. Grabisch et al. (2000) solved this minimization problem. Expressed in terms of the corresponding interaction parameters {β(λ) : λ ∈ Λ}, the optimal potential new parameter values are and the obtained minimum sum of squares is We use the latter to define the function d(λ ⋆ , τ, Λ, {θ(λ) : λ ∈ Λ}), used in (33) to define the distribution for what interaction λ ⋆ to remove. We simply set Combining the expression in (35) with the one-to-one relations in (9) and (10), one can find the potential new parameters {θ ⋆ (λ) : λ ∈ Λ ⋆ } in terms of the current parameters {θ(λ) : λ ∈ Λ}. In particular, we see that this relation is linear and we have a |Λ| × |Λ| matrix A so that where θ ⋆ = (θ ⋆ (λ) : λ ∈ Λ) T and θ = (θ(λ) : λ ∈ Λ) T are column vectors of the potential new and current parameter values, respectively. As the number of elements in θ ⋆ is one less than the number of elements in θ, we use β(λ ⋆ ) to obtain the one-to-one relation we need for a reversible jump proposal. The Jacobian determinant in the expression for the corresponding acceptance probability is clearly det(A), and in Appendix B we show that the absolute value of this determinant is always equal to one, i.e. | det(A)| = 1.

Proposing to add an inactive interaction to Λ
If it is decided that an inactive interaction should be added to Λ, the next step is to decide what interaction λ ⋆ ∈ Ω(τ 0 ) \ Λ to add. We do this in two steps, first we draw at random whether a first-order or a higher-order interaction should be added to Λ. If a first-order interaction should be added, we draw uniformly at random a node v ⋆ from τ 0 \ τ and set λ ⋆ = {v ⋆ }. Then τ ⋆ = τ ∪ λ ⋆ and Λ ⋆ = Λ ∪ {λ ⋆ }. If τ = τ 0 , so that no such v ⋆ exist, we simply propose an unchanged state. If a higher-order interaction should be added we need to ensure that Λ ∪ {λ ⋆ } is dense. We therefore first find and thereafter draw λ ⋆ uniformly at random from Λ a . Then τ ⋆ = τ and Λ ⋆ = Λ∪{λ ⋆ }. If it should happen that Λ a = ∅, we again simply propose an unchanged state.
Having decided τ ⋆ and Λ ⋆ , the next step is to generate the potential new parameter values {θ ⋆ (λ) : λ ∈ Λ ⋆ }. When doing this, one should remember that this adding a potential new interaction proposal must be one-to-one with the reverse removing an interaction proposal discussed in Section 5.2.1. Therefore, the proposal distribution for the potential new parameter values {θ ⋆ (λ) : λ ∈ Λ ⋆ } must conform with (35), and thereby also with (38). A natural way to achieve this is to draw a value β ⋆ (λ ⋆ ) from some distribution and define the potential new interaction parameters by the inverse transformation of (35), i.e.
The difference between what we now have to do and what is done in the Gibbs update is that in the Gibbs update the values ∆(λ) are sampled independently from a Gaussian distribution, whereas here these are implicitly defined by (40) together with the one-to-one relations (9) and (10). It is tempting to sample α = β ⋆ (λ ⋆ ) from the resulting full conditional, as this would give a high density for values of β ⋆ (λ ⋆ ) that corresponds to models with a high posterior probability.
As discussed in Section 5.1 for the Gibbs update, it is computationally feasible to sample from this full conditional by adaptive rejection sampling. However, the normalizing constant of this full conditional is not computationally available, and for computing the associated acceptance probability the normalizing constant of the distribution of β ⋆ (λ ⋆ ) must be available. To construct a proposal distribution for β ⋆ (λ ⋆ ) = α, we therefore instead first generate r (say) independent samples α 1 , . . . , α r from the full conditional for α, by adaptive rejection sampling, and thereafter draw α = β ⋆ (λ ⋆ ) from a Gaussian distribution with mean valueᾱ = 1 n r i=1 α i and variance s 2 α = 1 r−1 r i=1 (α i −ᾱ) 2 . Our proposal distribution for β ⋆ (λ ⋆ ) is thereby an approximation to its full conditional.
As this is a reversible jump proposal, the associated acceptance probability includes a Jacobian determinant. By construction the Jacobian determinant for this proposal is the inverse of the Jacobian determinant for the removing an interaction proposal discussed in Section 5.2.1. As we have | det(A)| = 1, we also get | det(A −1 )| = 1.

Simulation examples
In this section we investigate our prior and proposal distributions on two binary example scenes. Firstly, we consider a mortality map for liver and gallbladder cancers for white males from 1950 to 1959 in the eastern United States, compiled by Riggan et al. (1987). Using Markov random field models, this data set has previously been analyzed by Sherman et al. (2006), Liang (2010) and Austad and Tjelmeland (2017), see also Liang et al. (2011). Secondly, we consider a data set previously considered by Stien and Kolbjørnsen (2011). They also fitted a Markov mesh model to this data set, but with manually chosen neighborhood and interaction structures. In the following we first discuss some general aspects relevant for both the two examples and thereafter present details of each of the two examples in turn.
As also briefly discussed in Section 5, we reduce the boundary effects of the Markov mesh model by letting x include an unobserved boundary around the observed area. We choose the unobserved boundary large enough so that each of the observed nodes are at least 20 nodes away from the extended lattice boundary.
We let χ denote the set of nodes in the extended lattice and let x = (x v , v ∈ χ) be the corresponding collection of binary variables. We assume x to be distributed according to the Markov mesh model defined in Section 3, and for τ , Λ and {θ(λ) : λ ∈ Λ} we adopt the prior specified in Section 4. We let χ o ⊂ χ denote the set of nodes for which we have observed values. Thereby χ u = χ\ χ o is the set of unobserved nodes. Correspondingly, we let To simulate from this distribution, we adopt the updates discussed in Section 5 to update τ , Λ and {θ(λ), λ ∈ Λ} conditioned on x = (x o , x u ), and we use singlesite Gibbs updates for each unobserved node v ∈ χ u given τ , Λ, {θ(λ), λ ∈ Λ} and x χ\{v} . We define one iteration of the algorithm to include |χ u | single-site Gibbs updates for randomly chosen nodes in χ u followed by either one Gibbs update of the parameter values {θ(λ), λ ∈ Λ} as discussed in Section 5.1 or one update of the active interactions as discussed in Section 5.2. In each iteration we independently update the parameter values or the active interactions with probabilities 0.55 and 0.45 respectively.
The prior defined in Section 4 contains three hyper-parameters, the radius r which defines the set of possible neighbors, the probability p ⋆ in (25), and the parameter σ in (27). In both examples, we use r = 5 which gives the 34 possible neighbors shown in Figure 4. To get a prior where the probability for a Markov mesh model with higher-order interactions is reasonably high, we set the value of p ⋆ as high as 0.9, and to get an essentially uniform prior distribution for p(λ), we set σ = 100. The proposal distribution discussed in Section 5.2 has one algorithmic tuning parameter, ν, and based on simulation results in preliminary runs we set ν = 0.5.
In the following we present the example scene and discuss corresponding simulation results for each of our two examples. We start with the cancer mortality map compiled by Riggan et al. (1987).

Cancer mortality map
The cancer mortality map data are shown in Figure 6(a), where black (x v = 1) and white (x v = 0) pixels represent counties with high and low cancer mortality rates, respectively. The gray area around the observed map represents unobserved nodes which we included in the model to reduce the boundary effects of the Markov mesh model. Adopting the Markov mesh and prior models discussed in Sections 3 and 4, respectively, with the hyper-parameters defined above, we use the RJMCMC setup discussed above to explore the resulting posterior distribution. We run the Markov chain for 2 500 000 iterations, and study trace plots of different scalar quantities to evaluate the convergence and mixing properties To study the posterior distribution we first estimate, for each of the 34 apriori potential neighbors in τ 0 , the posterior probability for v ∈ τ 0 to be a neighbor.
To estimate this we simply use the fraction of simulated models where v is in the template sequential neighborhood τ . The result is shown in Figure 6 where the estimated probabilities are 0.049388 and 0.030353, respectively. From the data set shown in Figure 6(a), we see that the dependence between neighbor nodes seems to be quite weak, so the low number of simulated neighbors should come as no surprise.   Table with  We also estimate the a posteriori marginal distributions for the parameter values θ(·) corresponding to the four high probable interactions. Note that some of the interactions do not exist in some of the simulated models, but the θ(·) value is still well defined and can be computed as discussed in Section 2.2. Figure 8 depicts the histograms of the simulated parameter values θ(·). From the simulation we also estimate the posterior probability for each of the possible models. The two most probable models are shown in Figure 9. These two models have posterior probabilities as high as 0.475 and 0.381 while the remaining probability mass is spread out on a very large number of models.
Finally, we generate realizations from simulated Markov mesh models. Figure   10 contains realizations simulated from four randomly chosen models simulated in the Markov chain (after the specified burn-in). As in Figure 6(a), showing the observed data set, black and white nodes v represent x v = 1 and 0, respectively. Comparing the realizations with the data set in Figure 6(a), we can get a visual impression of to what degree the simulated models have captured the dependence structure in the data set. To study this also more quantitatively, we consider the 16 possible configurations in a 2 × 2 block of nodes. For each of these configurations, we find in a realization the fraction of such blocks that has the specified configuration. By repeating this for a large number of realizations we estimate the posterior distribution for the fraction of 2 × 2 blocks with a specified configuration in a realization. This distribution should be compared with the corresponding fraction in the observed data set. Figure 11 shows the estimated density for each of the 16 configurations. The corresponding fractions for the observed data set are marked by vertical dotted lines. Note that for most of these distributions the corresponding fractions for the observed data set are  centrally located in the distribution. The exceptions are (g) and partly (i) and (j), where the observed quantity is more in the tail of the distribution.

Sisim data set
In this example we reconsider a data set previously studied in Stien and Kolbjørnsen (2011). The scene, shown in Figure 12(a), is simulated by the sequential indicator simulation procedure (Journel, 1982;Deutsch and Journel, 1998) and it is a much used example scene in the geostatistical community. We name the data set "sisim". The sisim scene is represented on a 121 × 121 lattice. To reduce the boundary effects of the Markov mesh model we again include unobserved nodes around the observed area, shown as gray in Figure 12(a).
Again adopting the Markov mesh and prior models defined in Sections 3 and 4 and the hyper-parameters defined above, we use the RJMCMC setup discussed above to explore the resulting posterior distribution. For this data set each iteration of the algorithm requires more computation time than in the cancer mortality map data, so we run the Markov chain for only 1 250 000 iterations. To evaluate the convergence properties of the simulated Markov chain, we study trace plots of different scalar quantities in the same way as in Section 6.1. Figure 13 shows trace plots of the first 50 000 iterations for the number of interactions and for the logarithm of the posterior density. At first glance at these two and the other trace plots we have studied, we asserted that the simulated chain had converged at least within the first 30 000 − 40 000 iterations. As an extra precaution we discarded the first 250 000 iterations when estimating posterior properties.
As in Section 6.1, we estimate the posterior probability for v ∈ τ 0 to be in the template sequential neighborhood τ . The results are shown in Figure 12(b), where we use a gray scale to visualize the probabilities. There are five nodes whose estimated posterior probabilities are essentially equal to 1, and these are would not gain much extra information by including both of them in the template sequential neighborhood. Moreover, the prior prefers parsimonious models, which we obtain by not including too many nodes in the template sequential neighborhood.
Next, as for the cancer mortality map data set, we correspondingly estimate the posterior probabilities for each possible interaction to be included in the  Table 2 shows the top 20 a posteriori most likely interactions and corresponding estimated probabilities. We see that many interactions have high posterior probabilities.
We also estimate the a posteriori marginal distributions for the parameter values θ(·) corresponding to the top eight most likely interactions. Figure 14 depicts the histograms of the simulated parameter values θ(·). From the simulation we also estimate the posterior probability for each of the possible models. The most probable model is shown in Figure 15. This model has posterior probability equal to 0.13802. The remaining probability mass is spread out on a very large number of models.
As in the cancer mortality data set example, we also now generate realizations from the simulated Markov mesh models. Figure 16 contains realizations simulated from four randomly chosen models simulated in the Markov chain (after the specified burn-in). As in Figure 12(a), showing the observed data set, black and white nodes v represent x v = 1 and 0, respectively. Also now we estimate the distribution of values in a 2 × 2 block of nodes. Figure 17 shows the estimated density for each of the 16 configurations. The corresponding fractions for the ob-   In the cancer mortality data set example, essentially all of the posterior probability mass was concentrated in a few models. In the sisim data set example, the probability mass is spread out on a very large number of models. In particular, as also discussed above, the most probable model has a posterior probability giving a clear indication that the posterior model distribution have more than one mode. To find a second mode we limit the attention to the simulated models that was not included in the first cluster of models and repeat the process. Thus, we first find the a posteriori most probable model not included in the first model cluster. This model is shown in Figure 18. Then we find all neighbors of this model, all neighbors of the neighbors and so on. The estimated posterior probability in this second cluster of models is 0.146563. Thus, these two first clusters contain more than 95% of the simulated models, and we therefore choose not to search for a third cluster. Knowing that we have two important clusters or modes it is natural to reconsider the convergence and mixing properties of our Markov chain. Figure 19 shows a trace plot of the visited clusters for the subsampled models, where 1 and 2 on the y-axis represent the first and second clusters found, respectively, and 3 represent all remaining models. We then see that the second cluster is in fact visited only once, giving a clear indication of poor mixing. We should thereby not trust the estimated probabilities for the two clusters given above, but that the chain is first moving from the first cluster to the second and thereafter back again clearly shows that both of them have a significant posterior probability mass.

Closing remarks
In this article we propose a prior distribution for a binary Markov mesh model.
The specification of a Markov mesh model has three parts. First a sequential neighborhood is specified, next the parametric form of the conditional distributions is defined, and finally we assign values to the parameters. We formulate prior distributions for all these three parts. To favor parsimonious models, our prior in particular assigns positive prior probabilities for some interaction parameters to be exactly zero. A corresponding prior formulation has previously been proposed for Markov random fields (Arnesen and Tjelmeland, 2017). The advantage of using it for a Markov mesh model is that an explicit and easy to compute expression is available for the resulting posterior distribution, whereas the posterior based on a Markov random field will include the computationally intractable normalizing constant of the Markov random field.
To sample from the resulting posterior distribution when conditioning on an observed scene, we adopt the RJMCMC setup. We propose an algorithm based on the combination of two proposal distributions, a Gibbs proposal for the parameter values and a reversible jump proposal changing the sequential neighborhood and parametric form of the conditional distributions.
To explore the performance of the specified prior distribution and the corresponding RJMCMC posterior simulation algorithm, we consider two scenes.
The first is an observed cancer mortality map data with small spatial coupling between neighboring nodes. For this scene the RJMCMC algorithm converges quickly and has good mixing properties. Most of the posterior mass ends up in models with only two nodes in the sequential neighborhood. The second scene we tried is a frequently used scene in the geostatistical community. It has more spatial continuity than the first scene. The convergence of the RJMCMC algorithm becomes much slower when conditioning on this scene. In particular the posterior seems to have at least two modes and the mixing between the modes is slow. Our simulation results indicate that the a posteriori most likely model has six nodes in the sequential neighborhood and the conditional distributions has a parametric form with as much as twelve parameters. This shows that the specified prior is flexible in that the model complexity favored by the corresponding posterior adapts to the the complexity of the scene conditioned on.
In this article we have focused on binary Markov mesh models and thereby binary scenes. Our strategy for prior specification and posterior simulation, however, can easily be extended to a situation with more than two colors. The main challenge in this generalization does not lie in the specification of the prior, but is computational in that one should expect the convergence and mixing of a corresponding RJMCMC algorithm to be slower for a multi-color model. A direction for future research is therefore to improve the proposal distributions to obtain better convergence and mixing for the RJMCMC algorithm, both in the binary and multi-color cases. In particular we think a promising direction here is to define an MCMC algorithm where several Metropolis-Hastings proposals can be generated in parallel and where the proposals may have added and removed more than just one interaction relative to the current model.
As also touched on in Section 5.2.1, the transformation in (38) can be done in three steps. First θ is transformed into a vector β of the corresponding current interaction parameters {β(λ) : λ ∈ Λ}. This relation is given in (10) and is in particular linear so we can write Arranging also the vector β so that lower order interactions comes first, it is easy to see from (10) that A 1 is a lower triangular matrix with all diagonal elements equal to one. Thus det(A 1 ) = 1. The second step in the transformation is to use (35) to define a vector β ⋆ containing the set of potential new interaction parameters {β ⋆ (λ) : λ ∈ Λ ⋆ }. As the proposal is to remove an interaction, the number of elements in β ⋆ is one less than the number of elements in β. To obtain a one-to-one relation as required in the reversible jump setup, we include the current value β(λ ⋆ ) in a vector together with β ⋆ . We let β(λ ⋆ ) be the last element in the vector and we arrange also the vector β ⋆ so that lower order interaction parameters come first. As the relation in (35) is linear we can then where the elements of the square matrix A 2 is defined by (35). To find the determinant of A 2 , let r denote the number of elements in β before β(λ ⋆ ), so that element number r + 1 in β is β(λ ⋆ ). From (35) it then follows that A 2 has the block structure, where I r×r is the r ×r identity matrix, A 12 2 is an r ×(|Λ|−r) matrix, 0 (|Λ|−r)×r is a (|Λ| − r) × r matrix of only zeros, and A 22 2 is the (|Λ| − r) × (|Λ| − r) permutation matrix where the elements (i, i + 1) for i = 1, . . . , |Λ| − r − 1 and (|Λ| − r, 1) equals one and all other elements are zero. Thereby we have det(A 2 ) = det(I r×r )· det(A 22 2 ) = det(A 22 2 ), and as A 22 2 is a permutation matrix its determinant is plus or minus one. Thus, | det(A 2 )| = 1. The third step in the transformation from θ to θ ⋆ is to use (9) to transform the vector of potential new interaction parameters, β ⋆ , to a corresponding vector θ ⋆ of potential new parameter values. As the relation in (9) is also linear, we can write where the elements of the matrix A 3 is defined by (9). Recalling that we have arranged the elements in both θ ⋆ and β ⋆ so that parameters corresponding to lower order interactions come first, it is easy to see from (9) that A 3 is an upper triangular matrix with all diagonal elements equal to one. Thus det(A 3 ) = 1.