Prior specification for binary Markov mesh models
Håkon Tjelmeland Xin Luo
stat.ME stat.CO




Prior specification for binary
 Markov mesh models 
 
 
Xin Luo and Håkon Tjelmeland
Department of
    Mathematical Sciences, Norwegian University of Science and
    Technology













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.





Key words: Markov mesh model, prior construction, pseudo-Boolean functions, reversible jump MCMC, 
sequential neighborhood. 



   



§ 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
<cit.> and <cit.> for MRFs, and <cit.> and <cit.> 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 
<cit.>. In particular, <cit.> 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 
<cit.> 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 <cit.>. 
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 <cit.>. 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 hyper-prior, we adapt the general strategy used in 
<cit.> 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 <ref> we
introduce most of the notations we use for defining our Bayesian Markov mesh model, and in particular
discuss pseudo-Boolean functions. In Section <ref> we use this to formulate 
the Markov mesh model class. In Section <ref> we construct our prior
distribution, and in Section <ref> 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 <ref> we present two simulation examples and lastly we give 
some closing remarks in Section <ref>.




§ 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 ℤ^2, where ℤ={ 0,± 1,± 2,…} is the set of 
all integers, and we use v=(i,j)∈ℤ^2 also to denote a node in such an infinite
lattice. We use λ,λ^⋆⊆ℤ^2 to denote
arbitrary sets of nodes. To translate a node (i,j)∈ℤ^2 by an amount (k,l)∈χ,
we adopt the notation

    (i,j) ⊕ (k,l) = (i+k,j+l).

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 λ⊆ℤ^2 by the same amount (k,l)∈χ,
we write

    λ⊕ (k,l) = { (i,j) ⊕ (k,l): (i,j) ∈λ}.

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 so-called 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 (i,j).
We let ρ_v = ρ_(i,j)⊂χ denote the predecessor set of a node v=(i,j)∈χ, i.e.

    ρ_(i,j) = {(k,l)∈χ: nk+l < ni+j},

see the illustration in Figure <ref>(a).

We let ν_v=ν_(i,j)⊆ρ_(i,j) denote the sequential neighborhood for node v=(i,j)∈χ
as illustrated in Figure <ref>(b).
In Section <ref> 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 

    τ⊂ψ = { (i,j): i∈ℤ,j∈ℤ^-}∪{ (0,j): j∈ℤ^-},

where ℤ^- = { -1,-2,…} is the set of all negative integers. 
The sequential neighborhood of a node v∈χ is then 
defined as 

    ν_v = ( τ⊕ v) ∩χ.

As illustrated in Figure <ref>,
 
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 
x_v=x_(i,j)∈{0,1}. 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. 

    ξ(x) = { v∈χ: x_v=1}.

In particular, the set of nodes in the sequential neighborhood of node v that is on is then,
using (<ref>) and that ξ(x)⊆χ,

    ξ(x)∩ν_v = ξ(x)∩ (τ⊕ v).

In the next section, we define the class of pseudo-Boolean functions which we in Section <ref>
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 <cit.> and 
<cit.>, but in the present article we adopt the formulation in (<ref>) as this gives 
simpler expressions when formulating our Markov mesh model in Section <ref> and the 
corresponding prior distribution in Section <ref>.

<cit.> 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

    β(λ) = θ(λ) + ∑_λ^⋆⊂λ (-1)^|λ∖λ^⋆|θ(λ^⋆)

The one-to-one relation in (<ref>) and (<ref>) is known as Moebious inversion, see for example
<cit.>.

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 (<ref>)
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, 

    {{v}: v∈τ} = {λ∈Λ: |λ| = 1}.

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 <cit.>, the set of active interactions Λ can be visualized by a 
directed acyclic graph (DAG), where we have one vertex for each active interaction λ∈Λ and a 
vertex λ∈Λ is a 
child of another vertex λ^⋆∈Λ if and only if λ=λ^⋆∪{ v} for some 
v∈τ∖λ^⋆. Figure <ref>

shows such a DAG for Λ={∅, {(0,-1)}, {(-1,0)}, {(-1,-1)}, {(-1,1)}, {(0,-1),(-1,0)},
{(0,-1),(-1,1)}}, which is based on τ={(0,-1),(-1,-1),(-1,0),(-1,1)}. This τ
can be used to define the sequential neighborhoods for nodes in a rectangular lattice as 
discussed in Section <ref>. 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 ⊠.

As also discussed in <cit.>, one should note that a pseudo-Boolean function θ(·) that 
is represented on a dense set Λ⊆Ω(τ) can be uniquely specified by 
the values of {θ(λ):λ∈Λ}. The remaining values of the pseudo-Boolean function,
θ(λ),λ∈Ω(τ)∖Λ, are then given by (<ref>) and (<ref>)
and the restriction β(λ)=0 for λ∉Λ. Moreover, as the relations in 
(<ref>) and (<ref>) are linear, each θ(λ), λ∈Ω(τ)∖Λ is
a linear function of  {θ(λ):λ∈Λ}.





§ MARKOV MESH MODEL

In this section we formulate a homogeneous binary Markov mesh model <cit.> for a
rectangular m× n lattice. We adopt the notation introduced in Section <ref>, 
so in particular χ denotes the set of 
all nodes in the m× n lattice and x = (x_v,v∈χ) is the collection of 
the binary variables associated to χ. In a Markov mesh model the 
distribution of x is expressed as

    f(x) = ∏_v∈χ f(x_v|x_ρ_v),

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 

    f(x_v|x_ρ_v) = f(x_v|x_ν_v),

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 <ref>, we assume the sequential neighborhoods ν_v,v∈χ to be defined 
as translations of a template sequential neighborhood τ 
as specified in (<ref>) and (<ref>).
Using the result in (<ref>), the conditional distribution 
f(x_v|x_ν_v) can then be uniquely represented by a pseudo-Boolean function
θ_v(λ),λ⊆τ by the relation 

    f(x_v|x_ρ_v)=exp{ x_v·θ_v( ξ(x) ∩ (τ⊕ v))}/1 + exp{θ_v( ξ(x) ∩ (τ⊕ v))}.

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

    f(x) = ∏_v∈χexp{ x_v·θ( ξ(x) ∩ (τ⊕ v))}/1 + exp{θ( ξ(x) ∩ (τ⊕ v))}.

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 higher-order 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

    f(τ,Λ,{θ(λ): λ∈Λ}) = f(τ) f(Λ|τ) f({θ(λ):λ∈Λ}|τ,Λ),

where f(τ) is a prior for the template sequential neighborhood τ, f(Λ|τ) is a prior for the set of active 
interactions 
Λ when τ is given, and f({θ(λ):λ∈Λ}|Λ) 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 (<ref>). 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 <ref> we use this with r=5, see the illustration in Figure <ref>.


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

    f(τ) = f(n_τ) f(τ|n_τ).

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

    f(n_τ) = 1/n_τ + 1.

Moreover, having chosen τ to be uniform given n_τ=|τ|, we get

    f(τ|n_τ) = 1|τ_0|n_τ,

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) <cit.>, 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 <ref>, 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, for k=0,1,…,|τ| we define 

    Ω_k(τ) = {λ∈Ω(τ): |λ|=k}Λ_k = {λ∈Λ: |λ|=k}.

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 <ref> and in particular (<ref>). 
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 

    f(Λ|τ) = ∏_k=2^|τ| f(Λ_k |Λ_k-1).

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
λ∈Ω_k(τ) can then be active only if all 
k-1'th order interactions λ^⋆∈Ω_k-1(λ) are active. We let Π_k denote this set of 
possible active k'th order interactions, i.e. we must have

    Λ_k ⊆Π_k = {λ∈Ω_k(τ): λ^⋆∈Λ_k-1λ^⋆⊂λ}.

We assume each interaction λ∈Π_k to be active with some probability p_k, independently of each other, and get

    f(Λ_k|Λ_k-1)=p_k^|Λ_k|(1-p_k)^|Π_k|-|Λ_k|Λ_k⊆Π_k.

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

    p_k = 
    p^⋆    |Π_k| ≤ |Λ_k-1|,
    
    p^⋆·|Λ_k-1||Π_k|

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 (<ref>) and (<ref>),
we have that each 
θ(λ),λ∈Λ have a one-to-one correspondence with the conditional probability

    p(λ) = f(x_v=1|x_ρ_v) = exp{θ(λ)}/1+exp{θ(λ)}.

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 <cit.> 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 θ(λ),

    f(θ(λ)|τ,Λ) ∝e^θ(λ)(1+e^θ(λ))^2· e^-θ(λ)^2/2σ^2,

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 <ref>.

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 (<ref>) is 
required when updating Λ in a reversible jump Metropolis–Hastings algorithm targeting a corresponding posterior distribution, but since 
(<ref>) is a univariate distribution this normalizing constant can easily be found by numerical integration.
Letting c(σ) denote the normalizing constant of f(θ(λ)|τ,Λ) the complete expression for the prior for 
{θ(λ):λ∈Λ} is

    f({θ(λ):λ∈Λ}|τ,Λ)=
    ∏_λ∈Λ[ c(σ) ·e^θ(λ)(1+e^θ(λ))^2· e^-θ(λ)^2/2σ^2].


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 <ref>. We adopt the prior defined in Section <ref> and 
want to explore the resulting posterior distribution 

    f(τ,Λ,{θ(λ):λ∈Λ}|x) ∝ 
    f(τ,Λ,{θ(λ):λ∈Λ}) f(x|τ,Λ,{θ(λ):λ∈Λ}),

by a reversible jump Markov chain Monte Carlo algorithm (RJMCMC), see <cit.>. 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 <ref>, 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. 



 §.§ Gibbs update for the parameter values 
{θ(λ):λ∈Λ} 

Let τ, Λ and {θ(λ):λ∈Λ} be the current state. In this update, we keep τ and 
Λ unchanged and generate new parameter values {θ^⋆(λ):λ∈Λ}. 
To generate the new parameter values we first draw a random direction {Δ(λ):λ∈Λ}
by sampling Δ(λ) from a standard normal distribution, independently for each λ∈Λ. We then set

    θ^∗(λ) = θ(λ) + αΔ(λ),

where α∈ℝ is sampled from the full conditional 

    f(α|τ,Λ,{θ(λ) + αΔ(λ):λ∈Λ},x)    ∝
    f({θ(λ) + αΔ(λ):λ∈Λ}|τ,Λ) 
       · f(x|τ,Λ,{θ(λ)+αΔ(λ):λ∈Λ}).

As α is sampled from its full conditional, this is a Gibbs update and the Metropolis–Hastings acceptance 
probability is one. The full conditional (<ref>) for α is not of a standard form, but 
in Appendix <ref> we show that it is log-concave, so to generate samples from it we adopt
the adaptive rejection sampling algorithm of <cit.>. 








 §.§ 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. 





  §.§.§ Proposing to remove an active interaction from Λ

Having decided that an interaction should be removed, the next step is to decide what 
interaction λ^⋆∈Λ to remove. As the potential new Λ^⋆=Λ∖{λ^⋆}
should be dense, we first find the set of active interactions λ^⋆ that fulfill this requirement, 

    Λ_r={λ∈Λ∖{∅} : Λ∖{λ}}.

Thereafter we draw what interaction λ^⋆∈Λ_r to be removed, with probabilities

    q(λ^⋆) = exp{-ν d(λ^⋆,τ,Λ,{θ(λ):λ∈Λ})}/∑_λ∈Λ_rexp{ -ν d(λ,τ,Λ,{θ(λ):λ∈Λ})}λ^⋆∈Λ_r,

where ν≥ 0 is an algorithmic tuning parameter to be specified, and 
d(λ^⋆,τ,Λ,{θ(λ):λ∈Λ})
is a function that should measure the difference between the current pseudo-Boolean function defined by 
τ, Λ and {θ(λ):λ∈Λ} and the potential new pseudo-Boolean function 
defined by τ^⋆, Λ^⋆ and {θ^⋆(λ):λ∈Λ^⋆}. The precise
formula we use for d(λ^⋆,τ,Λ,{θ(λ):λ∈Λ}) we specify below,
after having specified how to set the potential new parameter values {θ^⋆(λ):λ∈Λ^⋆}.
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 (<ref>) and (<ref>).
Moreover, together with the restriction β(λ)=0 for λ∉Λ, this defines a
pseudo-Boolean function {θ(λ):λ∈Ω(τ_0)}. Correspondingly, there is a one-to-one relation 
between the potential new parameter values {θ^⋆(λ):λ∈Λ} and a set of potential 
new interaction parameters {β^⋆(λ):λ∈Λ^⋆}, and together with the restrictions
β^⋆(λ)=0 for λ∉Λ^⋆ this defines a potential new pseudo-Boolean function 
{θ^⋆(λ):λ∈Ω(τ_0)}. To get a high acceptance probability for the proposed change,
it is reasonable to choose the potential new parameter values {θ^⋆(λ):λ∈Λ^⋆} so that 
the difference between the two pseudo-Boolean functions {θ(λ):λ∈Ω(τ_0)} and 
{θ^⋆(λ):λ∈Ω(τ_0)} is small. One may consider the potential new pseudo-Boolean function 
{θ^⋆(λ):λ∈Ω(τ_0)} as an approximation to the current 
{θ(λ):λ∈Ω(τ_0)} and, adopting a minimum sum of squares criterion, minimize

    ({θ^⋆(λ):λ∈Λ^⋆}) = 
    ∑_λ∈Ω(τ_0)( θ^⋆(λ) - θ(λ))^2

with respect to {θ^⋆(λ):λ∈Ω(τ_0)}. <cit.> solved this  
minimization problem. Expressed in terms of the corresponding interaction parameters 
{β(λ): λ∈Λ}, the optimal potential new parameter values are 

    β^⋆(λ) = {[ β(λ) - 
    (-1/2)^|λ^⋆|-|λ|β(λ^⋆)                               ;                           β(λ)                                ].

and the obtained minimum sum of squares is

    min{({θ^⋆(λ):λ∈Λ})} = β(λ^⋆)/2^|λ^⋆|.

We use the latter to define the function d(λ^⋆,τ,Λ,{θ(λ):λ∈Λ}),
used in (<ref>) to define the distribution for what interaction λ^⋆ to remove. We simply set

    d(λ^⋆,τ,Λ,{θ(λ):λ∈Λ}) = β(λ^⋆)/2^|λ^⋆|.

Combining the  expression in (<ref>) with the one-to-one relations in (<ref>) and (<ref>), 
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 

    [[    θ^⋆; β(λ^⋆) ]] = A θ⇔θ = A^-1[ [    θ^⋆; β(λ^⋆) ]],

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 (A), and in Appendix <ref> we show that the absolute 
value of this determinant is always equal to one, i.e. |(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

    Λ_a = {λ∈Ω(τ_0)∖Λ: |λ| > 1 Λ∪{λ}}

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 <ref>. Therefore, the proposal distribution for the potential new
parameter values {θ^⋆(λ): λ∈Λ^⋆} must conform with (<ref>),
and thereby also with (<ref>). 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 
(<ref>), i.e. 

    β^⋆(λ) = {[ β(λ) + (-1/2)^|λ^⋆|-|λ|β^⋆(λ^⋆)                               ,;                            β(λ)                                 ].



It now just remains to specify from what distribution to sample β^⋆(λ^⋆).
The potential new parameter values {θ^⋆(λ):λ∈Λ^⋆} are linear functions
of β^⋆(λ^⋆), and by setting β^⋆(λ^⋆)=α it can be expressed as in 
(<ref>) for the Gibbs update. 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 (<ref>) together with the one-to-one
relations (<ref>) and (<ref>). 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 <ref> 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∑_i=1^rα_i and 
variance s_α^2 = 1/r-1∑_i=1^r (α_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 <ref>. As we have 
|(A)|=1, we also get |(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 <cit.>. 
Using Markov random field models, this data set has previously been analyzed by <cit.>, <cit.> and <cit.>, see 
also <cit.>. Secondly, we consider a data set 
previously considered by <cit.>. 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 <ref>, 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 <ref>, and for τ, Λ and 
{θ(λ) : λ∈Λ} we adopt the prior specified in Section <ref>. 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 x_o=(x_v, v∈χ_o) be the observed values and 
x_u=(x_v,v∈χ_u) the unobserved values. The posterior distribution of interest is thereby 
f(τ,Λ,{θ(λ),λ∈Λ}|x_o). To simplify the posterior simulation, we include x_u as 
auxiliary variables and adopt the reversible jump Metropolis–Hastings algorithm to simulate from 

    f(τ,Λ,{θ(λ),λ∈Λ},x_u|x_o) ∝
    f(τ,Λ,{θ(λ),λ∈Λ}) f(x_o,x_u|τ,Λ,{θ(λ),λ∈Λ}).

To simulate from this distribution, we adopt the updates discussed in Section <ref> to update
τ, Λ and {θ(λ),λ∈Λ} conditioned on x=(x_o,x_u), and we use single-site 
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 <ref> or one update of the active 
interactions as discussed in Section <ref>. 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 <ref> contains three hyper-parameters, the radius r which defines the set of 
possible neighbors, the probability p^⋆ in (<ref>), and the parameter σ in (<ref>).
In both examples, we use r=5 which gives the 34 possible neighbors shown in Figure <ref>.
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 <ref> 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 <cit.>. 



 §.§ Cancer mortality map

The cancer mortality map data are shown in Figure <ref>, 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 <ref> and <ref>, 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 of the simulated Markov chain.
Figure <ref> 

shows trace plots of the first 25 000 iterations for the
number of interactions and for the logarithm of the posterior density.
From these two and the other trace plots we have studied, we conclude that the simulated chain has converged
at least within the first 10 000-15 000 iterations. As an extra precaution we discard the first 
25 000 iterations when estimating posterior 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 <ref>, where we use a gray scale to visualize the probabilities.
Nodes (0,-1) and (-1,0) have high estimated posterior probabilities, equal to 0.999819 and 0.990577, respectively. 
The third and fourth most probable neighbor nodes are (-1,-1) and (-1,2), where the estimated 
probabilities are 0.049388 and 0.030353, respectively. From the data set shown in Figure <ref>, 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.

Next we correspondingly estimate the posterior probabilities for each possible interaction to be included in the model. 
Table <ref> shows the top 10 a posteriori most likely interactions and the corresponding estimated probabilities. 

We see that the first four interactions have high posterior probabilities while the others have low probabilities. 
In addition, the four most likely interactions only include the high probability neighbor nodes (0,-1) and (-1,0).

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 <ref>. 
Figure <ref> 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 <ref>. 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 <ref> contains realizations 
simulated from four randomly chosen models simulated in the Markov chain (after the specified burn-in).

As in Figure <ref>, 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 <ref>, 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 <ref> 

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 
<cit.>. The scene, shown in Figure <ref>, 

is simulated by the sequential indicator
simulation procedure <cit.> 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 
<ref>. 


Again adopting the Markov mesh and prior models defined in Sections <ref> and <ref> 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 <ref>.
Figure <ref> 

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 <ref>, we estimate the posterior probability for v∈τ_0 to be in the template
sequential neighborhood τ. The results are shown in Figure <ref>, 
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 (0,-1), (-1,0), (-1,2), (0,-3) and (-1,4). Four more nodes have 
estimated posterior probabilities higher than 0.1. These are (-2,3), (-3,-1), (-2,-3) and (-1,3) with estimated probabilities 0.444608, 0.425779, 0.323181 and 0.182879, respectively. It is interesting to note the spatial locations of the high probability nodes. At least for a part of the area
every second node is chosen as a neighbor with high probability. To understand this effect, we must remember that the values of
two nodes that are lying next to each other are highly correlated, so one 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 model.
Table <ref> 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 <ref> 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 <ref>. 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 <ref> contains 
realizations simulated from four randomly chosen models simulated in the Markov chain (after the specified burn-in).

As in Figure <ref>, 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 <ref> 

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 (c), and partly (e), (f), (i) and (j), where the observed quantities are more in the tail of the distribution.


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 estimated to be as 
low as 0.13802. Using the 
simulated models to understand the posterior model distribution is then more difficult. As a first step in describing
the posterior model distribution, our focus here is on whether it has one or several modes. To do this we 
first need to define what we should mean by a mode in this complicated model space. We start by defining two models to be 
neighbors if one of them can be obtained from the other by including one extra interaction. Thus, our proposal
distribution in Section <ref>, proposing to change the set of active interactions, is 
always generating a potential new model that is a neighbor of the current model. To explore whether we have several modes
in our posterior distribution, we first subsample the simulated Markov chain, keeping a realization every
50 iterations after the burn-in period. This leave us with 
20 000 realizations.  From these we first find the most frequent model, 
visualized in Figure <ref>, and then all neighbor models to this most probable model, all 
neighbor models to the neighbors, and so on until the process stops. The sum of the estimated 
posterior probabilities of the models in the resulting cluster of models is 0.80755,
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 <ref>. 

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 <ref>  

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
<cit.>. 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.



jasa







§ LOG-CONCAVITY OF THE FULL CONDITIONAL FOR Α

In this appendix we prove that the full conditional f(α|τ,Λ,{θ(λ)+αΔ(λ):λ∈Λ},x)
defined in (<ref>) is log-concave, so that we can use adaptive rejection sampling to generate samples from it. Defining
g(α) = ln[ f(α|τ,Λ,{θ(λ)+αΔ(λ):λ∈Λ},x)]
and using (<ref>) we have

    g(α)    = ln[  f({θ(λ) + αΔ(λ):λ∈Λ}|τ,Λ)]
       + ln[ f(x|τ,Λ,{θ(λ)+αΔ(λ):λ∈Λ})]  + C,

where C is the logarithm of the normalizing constant in (<ref>). Inserting expressions for 
the prior and likelihood in (<ref>) and (<ref>), respectively, we get

    g(α)    =  ∑_λ∈Λ[ c(λ) +θ(λ)+αΔ(λ) 
    -2ln(1+e^θ(λ)+αΔ(λ)) - 
    (θ(λ)+αΔ(λ))^2/2σ^2] 
       + ∑_v∈χ[ x_v (θ(ξ(x)∩(τ⊕ v)) + αΔ(ξ(x)∩(τ⊕ v))) . 
       - .ln(1+e^θ(ξ(x)∩(τ⊕ v)) + αΔ(ξ(x)∩(τ⊕ v)))]
    + C.

Grouping terms of the same functional form, we get

    g(α)    = C_0 + C_1 α - 1/2σ^2∑_λ∈Λ (θ(λ)+αΔ(λ))^2 
    - 2∑_λ∈Λln(1+e^θ(λ)+αΔ(λ))
       - ∑_v∈χln(1+e^θ(ξ(x)∩(τ⊕ v)) + αΔ(ξ(x)∩(τ⊕ v))),

where 

    C_0 = C + ∑_λ∈Λθ(λ) +∑_v∈χx_vθ(ξ(x)∩ (τ⊕ v)) 
    C_1 = ∑_λ∈ΛΔ(λ) + ∑_v∈χΔ(ξ(x)∩ (τ⊕ v))

are constants as a function of α. The second derivative of the constant and linear terms in 
(<ref>) are of course zero. Since the coefficients of the quadratic terms are all negative, the second 
derivative of all of these are less or equal to zero, and unless Δ(λ) equals zero for 
all λ∈Λ the second derivative of the sum of these terms is even strictly less than zero. 
The remaining terms in (<ref>)
all have the same functional form 
as a function of α, namely

    h(α) = - a ln(1+e^b+cα),

a term in the sum over λ∈Λ has a=2, b=θ(λ) and c=Δ(λ), whereas a term in the 
sum over v∈χ has a=1, b=θ(ξ(x)∩ (τ⊕ v)) and c=Δ(ξ(x)∩ (τ⊕ v)). To prove 
that the second derivative of all of these terms are negative, it is thereby sufficient to show that h^''(α) < 0
for all a>0 and α,b,c∈ℝ. Simple differentiation gives

    h^''(α) = - ac^2 e^b+cα/(1+e^b+cα)^2.

Thus, h^''(α)<0 for all a>0 and α,b,c∈ℝ, and thereby g(α) is concave and the full conditional 
f(α|τ,Λ,{θ(λ)+αΔ(λ):λ∈Λ},x) is log-concave.




§ JACOBIAN DETERMINANT FOR THE PROPOSAL IN 
SECTION <REF>


The Jacobi determinant for our removing an active interaction from Λ proposal is (A), where A is defined by
(<ref>). The exact form of the matrix A depends on how we define the vectors θ and θ^⋆ used in (<ref>).
The vector θ should contain the set of current parameters {θ(λ):λ∈Λ}, but so far we 
have not specified what order to use when arranging this set of parameters into the vector θ. Correspondingly,
we have not specified what order to use when arranging the set of potential new parameters 
{θ^⋆(λ):λ∈Λ^⋆} into the vector θ^⋆. However, even if the elements of 
A depends on how we construct θ and θ^⋆, the absolute value of the determinant of A is the same 
for all arrangements of θ and θ^⋆.
To find (A) we arrange the vector θ so that parameters corresponding to lower order interactions comes first. 
The first element of the vector θ is thereby θ(∅), thereafter follows parameters corresponding to 
the first order interactions {θ(λ):λ∈Λ, |λ|=1} (in an arbitrary order), then all parameters
corresponding to second order interactions {θ(λ):λ∈Λ, |λ|=2} (again in an arbitrary order), 
and so on. We arrange θ^⋆ correspondingly, parameters corresponding to lower order interactions comes first.

As also touched on in Section <ref>, the transformation in (<ref>) can be done in three steps.
First θ is transformed into a vector β of the corresponding current interaction parameters 
{β(λ):λ∈Λ}. This relation is given in (<ref>) and is in particular linear so we can 
write

    β = A_1θ.

Arranging also the vector β so that lower order interactions comes first, it is easy to see from 
(<ref>) that A_1 is a lower triangular matrix with all diagonal elements equal to one. Thus (A_1)=1.
The second step in the transformation is to use (<ref>) 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 (<ref>) is linear we can then write 

    [[    β^⋆; β(λ^⋆) ]] = A_2β,

where the elements of the square matrix A_2 is defined by (<ref>). 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 (<ref>) it then follows that A_2 has the block structure,

    A_2 = [ [       I_r× r       A_2^12; 0_(|Λ|-r)× r       A_2^22 ]],

where I_r× r is the r× r identity matrix, A_2^12 is an r× (|Λ|-r) matrix, 0_(|Λ|-r)× r
is a (|Λ|-r)× r matrix of only zeros, and A_2^22 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 (A_2) = (I_r× r) ·(A_2^22) = (A_2^22), and as A_2^22 is a
permutation matrix its determinant is plus or minus one. Thus, |(A_2)|=1. The third step in 
the transformation from θ to θ^⋆ is to use (<ref>) to transform the vector of potential new 
interaction parameters, β^⋆, to a corresponding vector θ^⋆ of potential new parameter values. As the relation in 
(<ref>) is also linear, we can write

    [[    θ^⋆; β(λ^⋆) ]] = A_3
    [[    β^⋆; β(λ^⋆) ]],

where the elements of the matrix A_3 is defined by (<ref>). 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 (<ref>) that A_3 is an upper triangular
matrix with all diagonal elements equal to one. Thus (A_3)=1. Setting the 
three steps in the transformation together we have A=A_1 A_2 A_3 and thereby |(A)| = |(A_1)| 
· |(A_2)|· |(A_3)| = 1.





