Julien Stoehr

A review on statistical inference methods for discrete Markov random fields

Keywords: statistics, Markov random fields, parameter estimation, model selection

Developing satisfactory methodology for the analysis of Markov random field is a very challenging task. Indeed, due to the Markovian dependence structure, the normalizing constant of the fields cannot be computed using standard analytical or numerical methods. This forms a central issue for any statistical approach as the likelihood is an integral part of the procedure. Furthermore, such unobserved fields cannot be integrated out and the likelihood evaluation becomes a doubly intractable problem. This report gives an overview of some of the methods used in the literature to analyse such observed or unobserved random fields.

Introduction

The problem of developing satisfactory methodology for the analysis of spatially correlated data has been of a constant interest for more than half a century now. Constructing a joint probability distribution to describe the global properties of such data is somewhat complicated but the difficulty can be bypassed by specifying the local characteristics via conditional probability instead. This proposition has become feasible with the introduction of Markov random fields (or Gibbs distribution) as a family of flexible parametric models for spatial data (the Hammersley-Clifford theorem, [START_REF] Besag | Spatial Interaction and the Statistical Analysis of Lattice Systems (with Discussion)[END_REF]. Markov random fields are spatial processes related to lattice structure, the conditional probability at each nodes of the lattice being dependent only upon its neighbours, that is useful in a wide range of applications. In particular, hidden Markov random fields offer an appropriate representation for practical settings where the true state is unknown. The general framework can be described as an observed data y which is a noisy or incomplete version of an unobserved discrete latent process x.

Gibbs distributions originally appears in statistical physics to describe equilibrium state of a physical systems which consists of a very large number of interacting particles such as ferromagnet ideal gases [START_REF] Lanford | Observables at infinity and states with short range correlations in statistical mechanics[END_REF]. But they have since been useful in many other modelling areas, surged by the development in the statistical community since the 1970's. Indeed, they have appeared as convenient statistical model to analyse different types of spatially correlated data. Notable examples are the autologistic model [START_REF] Besag | Spatial Interaction and the Statistical Analysis of Lattice Systems (with Discussion)[END_REF] and its extension the Potts model. Shaped by the development of [START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF] and [START_REF] Besag | On the Statistical Analysis of Dirty Pictures[END_REF], -see for example [START_REF] Alfò | A finite mixture model for image segmentation[END_REF] and [START_REF] Moores | Segmentation of cone-beam CT using a hidden Markov random field with informative priors[END_REF] who performed image segmentation with the help of this modelling -and also in other applications including disease mapping (e.g., [START_REF] Green | Hidden Markov Models and Disease Mapping[END_REF] and genetic analysis (e.g., [START_REF] François | Bayesian Clustering Using Hidden Markov Random Fields in Spatial Population Genetics[END_REF][START_REF] Friel | Bayesian Inference in Hidden Markov Random Fields for Binary Data Defined on Large Lattices[END_REF] to name a few. The exponential random graph model or p * model [START_REF] Wasserman | Logit models and logistic regressions for social networks: I. An introduction to Markov graphs and p *[END_REF] is another prominent example [START_REF] Frank | Markov graphs[END_REF] and arguably the most popular statistical model for social network analysis (e.g., [START_REF] Robins | An introduction to exponential random graph (p*) models for social networks[END_REF].

Interests in these models is not so much about Markov laws that may govern data but rather the flexible and stabilizing properties they offer in modelling. Whilst the Gibbs-Markov equivalence provides an explicit form of the joint distribution and thus a global description of the model, this is marred by a considerable computational curse. Conditional probabilities can be easily computed, but the joint and the marginal distribution are meanwhile unavailable since the normalising constant is of combinatorial complexity and generally can not be evaluated with standard analytical or numerical methods. This forms a central issue in statistical analysis as the computation of the likelihood is an integral part of the procedure for both parameter inference and model selection. Remark the exception of small latices on which we can apply the recursive algorithm of [START_REF] Reeves | Efficient recursions for general factorisable models[END_REF], [START_REF] Friel | Recursive computing and simulation-free inference for general factorizable models[END_REF] and obtain an exact computation of the normalizing constant. However, the complexity in time of the aforementioned algorithm grows exponentially and is thus helpless on large lattices. Many deterministic or stochastic approximations have been proposed for circumventing this difficulty and developing methods that are computationally efficient and accurate is still an area of active research. Solutions to deal with the intractable likelihood are of two kinds. On one hand, one can rely on pseudo-model as surrogates for the likelihood. Such solutions typically stems from composite likelihood or variational approaches. On the other hand Monte Carlo methods have played a major role to estimate the intractable likelihood in both frequentist and Bayesian paradigm.

The present survey paper cares about the problem of carrying out statistical inference (mostly in a Bayesian framework) for Markov random fields. When dealing with hidden random fields, the focus is solely on hidden data represented by discrete models such as the Ising or the Potts models. Both are widely used examples and representative of the general level of difficulty. Aims may be to infer on parameters of the model or on the latent state x. The paper is organised as follows: it begins by introducing the existence of Markov random fields with some specific examples (Section 2). The difficulties inherent to the analysis of such a stochastic model are especially pointed out in Section 3. The special case of small regular lattices is described in Section 4. As befits a survey paper, Section 5 focuses on solutions based on pseudo-models while Section 6 is dedicated to a state of the art concerning sampling method.

Markov random field and Gibbs distribution 2.1 Gibbs-Markov equivalence

A discrete random field X is a collection of random variables X i indexed by a finite set S = {1, . . . , n}, whose elements are called sites, and taking values in a finite state space X := {0, . . . , K -1}. For a given subset A ⊂ S , X A and x A respectively define the random process on A, i.e., , {X i , i ∈ A}, and a realization of induces a topology on S : by definition, sites i and j are adjacent or neighbour if and only if i and j are linked by an edge in G . A random field X is a Markov random field with respect to G , if for all configuration x and for all sites i it satisfies the following Markov property

c 0 = , c 1 = , c 2 = C 4 = (c) c 3 = , c 4 = , c 5 = c 6 = , c 7 = , c 8 = c 9 = C 8 = C 4 ∪ (d)
P (X i = x i | X -i = x -i ) = P X i = x i X N (i) = x N (i) , (2.1)
where N (i) denotes the set of all the adjacent sites to i in G . It is worth noting that any random field is a Markov random field with respect to the trivial topology, that is the cliques of G are either the empty set or the entire set of sites S . Recall a clique c in an undirected graph G is any single vertex or a subset of vertices such that every two vertices in c are connected by an edge in G . As an example, when modelling a digital image, the lattice is interpreted as a regular 2D-grid of pixels and the random variables states as shades of grey or colors. Two widely used adjacency structures are the graph G 4 (first order lattice), respectively G 8 (second order lattice), for which the neighbourhood of a site is composed of the four, respectively eight, closest sites on a two-dimensional regular lattice, except on the boundaries of the lattice, see Figure 1.

The difficulty with the Markov formulation is that one defines a set of conditional distributions which does not guarantee the existence of a joint distribution. The Hammersley-Clifford theorem states that if the distribution of a Markov random field with respect to a graph G is positive for all configuration x then it admits a Gibbs representation for the same topology (see e.g., [START_REF] Grimmett | A theorem about random fields[END_REF], [START_REF] Besag | Spatial Interaction and the Statistical Analysis of Lattice Systems (with Discussion)[END_REF] and for a historical perspective [START_REF] Clifford | Markov random fields in statistics[END_REF]), namely a density function on X parametrised by θ ∈ Θ ⊂ R d and given with respect to the counting measure by

f (x | θ, G ) = 1 Z (θ, G ) exp {-V θ,G (x)} := 1 Z (θ, G ) q(x | θ, G ), (2.2)
where V θ,G denotes the energy function which can be written as a sum over the set C of all cliques of the graph, namely V θ,G (x) = c∈C V θ,c (x) for all configuration x ∈ X . The inherent difficulty of all these models arises from the intractable normalizing constant, sometimes called the partition function, defined by

Z(θ, G ) = x∈X exp {-V θ,G (x)} .
The latter is a summation over the numerous possible realisations of the random field X, that is of combinatorial complexity and cannot be computed directly (except for small grids and small number of states K). For binary variables X i , the number of possible configurations reaches 2 n .

Autologistic model and related distributions

The formulation in terms of potential allows the local dependency of the Markov field to be specified and leads to a class of flexible parametric models for spatial data. In most cases, cliques of size one (singleton c 0 ) and two (doubletons C = {c 1 , c 2 , c 3 , c 4 }) are assumed to be satisfactory to model the spatial dependency and potential functions related to larger cliques are set to zero. We present below popular examples of such pairwise models broadly used in the literature.

Autologistic model-Ising model

The autologistic model first proposed by [START_REF] Besag | Nearest-neighbour Systems and the Auto-Logistic Model for Binary Data[END_REF] is a pairwise-interaction Markov random field for binary (zero-one) spatial process. Denote θ = (α, β). The joint distribution is given by the following energy function

V θ,G (x) = -α n i=1 x i -β i G ∼j x i x j , (2.3)
where the above sum

i G ∼j
ranges the set of edges of the graph G . The full-conditional probability thus writes like a logistic regression where the explanatory variables are the neighbours and themselves observations. The parameter α controls the level of 0 -1 whereas the parameters {β} model the dependency between two neighbouring sites i and j. One usually prefers to consider variables taking values in {-1, 1} instead of {0, 1} since it offers a more parsimonious parametrisation and avoids non-invariance issues when one switches states 0 and 1 as mentioned by [START_REF] Pettitt | Efficient calculation of the normalizing constant of the autologistic and related models on the cylinder and lattice[END_REF]. A well known example is the general Ising model of ferromagnetism [START_REF] Ising | Beitrag zur theorie des ferromagnetismus[END_REF] that consists of discrete variables representing spins of atoms. The Gibbs distribution (2.3) is sometimes referred to as the Boltzmann distribution in statistical physics. The potential on singletons describes local contributions from external fields to the total energy. Spins most likely line up in the same direction of α, that is, in the positive, respectively negative, direction if α > 0, respectively α < 0. Putting differently α adjusts non-equal abundances of the two state values. The parameters {β} represent the interaction strength between neighbours i and j. When β > 0 the interaction is called ferromagnetic and adjacent spins tend to be aligned, that is neighbouring sites with same sign have higher probability. When β < 0 the interaction is called anti-ferromagnetic and adjacent spins tend to have opposite signs.

Potts model

The Potts model [START_REF] Potts | Some generalized order-disorder transformations[END_REF] originally appears in statistical mechanics to model interacting spins but has been used in other modelling area since then. It is a pairwise Markov random field that extends the Ising model to K possibles states. The model sets a probability distribution on parametrized by θ = (α 0 , . . . , α K-1 , β), namely

V θ,G (x) = - n i=1 K-1 k=0 α k 1{x i = k} -β i G ∼j 1{x i = x j }, (2.4)
where 1{A} is the indicator function equal to 1 if A is true and 0 otherwise. Note that a potential function can be defined up to an additive constant. To ensure that potential functions on singletons are uniquely determined, one usually imposes the constraint

K-1 k=0 α k = 0.

Phase transition

One major peculiarity of Markov random field in the absence of an external field (i.e., α = 0) is a symmetry breaking for large values of parameter β due to a discontinuity of the partition function when the number of sites n tends to infinity. When the parameter β is zero, the random field is a system of independent uniform variables and all configurations are equally distributed. Increasing β favours the variable X i to be equal to the dominant state among its neighbours and leads to patches of like-valued variables in the graph, such that once β tends to infinity values x i are all equal. The distribution thus becomes multi-modal. In physics this is known as phase transition. This transition phenomenon has been widely study in both physics and probability, see for example [START_REF] Georgii | Gibbs Measures and Phase Transitions[END_REF] for further details. The two dimensional Ising model is known to have a phase transition at a critical value β c . [START_REF] Onsager | Crystal Statistics. I. A Two-Dimensional Model with an Order-Disorder Transition[END_REF] obtained an exact value of β c for the Ising model on the first order square lattice, namely

β c = 1 2 log 1 + √ 2 ≈ 0.44.
The latter extends to a K-states Potts model on the first order lattice

β c = log 1 + √ K ,
see for instance [START_REF] Matveev | Complex-temperature singularities in Potts models on the square lattice[END_REF] for specific results to Potts model on the square lattice and [START_REF] Wu | The Potts model[END_REF] for a broader overview.

The transition is more rapid than the number of neighbours increases. To illustrate this point, Figure 2 gives the average proportion of homogeneous pairs of neighbours, and the corresponding variance, for 2-states Potts model on the first and second order lattices of size 100 × 100. Indeed, phase transition corresponds to discontinuity at

β c of β → lim n→∞ 1 n ∇ log Z (β, G ). Since ∇ 2 log Z(β, G ) = Var {S(X)}, where S(X) = i G ∼j 1{X i = X j }
is the number of homogeneous pairs of a Potts random field X, the discontinuity condition can thus be written as

lim β→βc lim n→∞ Var {S(X)} = ∞.
Mention this is all theoritical asymptotic considerations and the discontinuity does not show itself on finite lattice realizations but the variance becomes increasingly sharper as the size grows.

Hidden Gibbs random field

Hidden Markov random fields has encountered a large interest over the past decade. It offer an appropriate representation for practical settings where the true state is unknown and observed indirectly through another field; this permits the modelling of noise that may happen upon many concrete situations: image analysis, (e.g., [START_REF] Besag | On the Statistical Analysis of Dirty Pictures[END_REF], Forbes and Peyrard, 2003[START_REF] Alfò | A finite mixture model for image segmentation[END_REF][START_REF] Moores | Segmentation of cone-beam CT using a hidden Markov random field with informative priors[END_REF], disease mapping (e.g., [START_REF] Green | Hidden Markov Models and Disease Mapping[END_REF], genetic analysis [START_REF] François | Bayesian Clustering Using Hidden Markov Random Fields in Spatial Population Genetics[END_REF], to name but a few. The unobserved data is modelled as a discrete Markov random field X associated to an energy function V θ,G , as defined in (2.2). Given the realization x of the latent, the observation y is a family of random variables indexed by the set of sites S , and taking values in a set Y , i.e., y = (y i ; i ∈ S ), and are commonly assumed as independent draws that form a noisy version of the hidden field. Consequently, we set the conditional distribution of Y knowing X = x, also called emission distribution, as the product

g φ (y | x) = i∈S g i (y i | x i , φ) ,
where g i is the marginal noise distribution parametrized by φ, that is given for any site i. Those marginal distributions are for instance discrete distributions [START_REF] Everitt | Bayesian Parameter Estimation for Latent Markov Random Fields and Social Networks[END_REF], Gaussian (e.g., [START_REF] Besag | Bayesian image restoration, with two applications in spatial statistics[END_REF], Forbes and Peyrard, 2003[START_REF] Cucala | Bayesian Inference on a Mixture Model With Spatial Dependence[END_REF] or Poisson distributions (e.g., [START_REF] Besag | Bayesian image restoration, with two applications in spatial statistics[END_REF]. Model of noise that takes into account information of the nearest neighbours have also been explored [START_REF] Besag | On the Statistical Analysis of Dirty Pictures[END_REF].

Assuming that all the marginal distributions g i are positive, one may write the joint distribution of (X, Y), also called the complete likelihood, as

P (x, y | θ, φ, G ) = 1 Z (θ, G ) exp -V θ,G (x) + i∈S log g i (y i | x i , φ) .
The conditional field X given Y = y is thus a Markov random field and the noise can be interpreted as a non homogeneous external potential on singleton which is a bond to the unobserved data.

Statistical analysis issues

The intractable normalising constant Z(θ, G ) forms a central issue for both parameter and model selection problems as the likelihood is an integral part of the statistical procedure. Below, we introduce some of the classical problems of the literature Maximum likelihood estimator Under the statistical model f (x | θ, G ), computing the maximum likelihood estimator, namely θMLE = arg max

θ log f (x | θ, G ) , (3.1)
is challenging. Indeed, closed-form gradients are typically out of reach. Furthermore, one cannot rely on differentiation techniques such as automatic differentiation (e.g., Neidinger, 2010) since point-wise estimation is impossible due to the intractability issue.

Computations of posterior distributions

Consider a Bayesian posterior distribution expressed as

π (θ | x) ∝ f (x | θ, G ) π(θ), (3.2)
where f (x | θ, G ) denoted the likelihood of the observed data x ∈ X and π(θ) denotes a prior density on the parameter space Θ with respect to a reference measure (often the Lebesgue measure of the Euclidean space). Here we are concerned with the situation where the unnormalised posterior distribution, the right-hand-side of (3.2) is intractable. This complication results in what is often termed a doubly-intractable posterior distribution, since the posterior distribution itself is normalised by the evidence (or marginal likelihood) which is typically also intractable.

Bayesian model choice

, . . . , ψ M ) given Y = y, π(m | y) = e(y | m)π(m) M m =1 e(y | m )π(m )
, where e(y | m) denotes the evidence (or integrated likelihood) of model m defined as

e(y | m) = Ψm f m (y | ψ m )π m (ψ m )dψ m . (3.3)
When the goal of the Bayesian analysis is the selection of the model that best fits the observed data y, it is performed through the maximum a posteriori (MAP) defined by

m MAP (y) = arg max m π(m | y). (3.4)
The standard approach to compare one model against another is based on the Bayes factor [START_REF] Kass | Bayes factors[END_REF] that involves the ratio of the evidence (3.3) of each model. Under the assumption of model being equally likely a priori, the MAP rule (3.4) is equivalent to choose the model with the largest evidence (3.3) and in place of a fully Bayesian approach, model choice criterion can be used. However, the evidence can usually not be computed with standard procedure because of a high-dimensional integral. Thus much of the research in model selection area focuses on evaluating it by numerical methods.

Selecting the model among a collection of Markov random fields is a daunting task as none of (3.3) and (3.4) are analytically available. The model selection problem for hidden Markov random fields is even more complicated and could be termed as a triply-intractable problem. Indeed in addition to the integral on Ψ m which is typically intractable, the stochastic model for Y is based on the latent process X in X , that is

f m (y | ψ m = (θ m , φ m )) = X g φ m (y | x)f (x | θ m , G m )µ(dx), (3.5)
with µ the counting measure (discrete case). Both the integral and the Gibbs distribution are intractable and consequently so is the posterior distribution.

Recursive algorithm for small Markov random field

When the Markov random field is defined on a small enough lattices, it is possible to answer the difficulty of computing the normalising constant by relying on generalised recursions for general factorisable models [START_REF] Reeves | Efficient recursions for general factorisable models[END_REF]. This method is based on an algebraic simplification due to the reduction in dependence arising from the Markov property. It applies to unnormalized likelihoods that can be expressed as a product of factors, each of which is dependent on only a subset of the lattice sites. More specifically the unnormalized version of a Gibbs distribution can be write as

q (x | θ, G ) = n-r i=1 q i (x i:i+r | θ, G ) ,
where each factor q i depends on a subset x i:r = {x i , . . . , x i+r } of x, where r is defined to be the lag of the model.

As a result of this factorisation, the summation for the normalizing constant can be represented as

Z (θ, G ) = xn-r:n q n-r (x n-r:n | θ, G ) • • • x1:1+r q 1 (x 1:1+r | θ, G ) .
The latter can be computed much more efficiently than the straightforward summation over the K n possible lattice realisations using the following steps

Z 1 (x 2:1+r ) = x1 q 1 (x 1:1+r ) , Z i (x i+1:i+r ) = xi q i (x i:i+r ) Z i-1 (x i:i+r-1 ) , for all i ∈ {2, . . . , n -r}, Z (ψ, G ) = xn-r+1:n Z n-r (x n-r+1:n ) .
The complexity of the troublesome summation is significantly cut down since the forward algorithm solely relies on K r possible configurations. Consider a rectangular lattice h × w = n, where h stands for the height and w for the width of the lattice, with a first order neighbourhood system G 4 (see Figure 1.(a)). The minimum lag representation for a pairwise model define on such a lattice occurs for r given by the smaller of the number of rows or columns in the lattice.

Without the loss of generality, assume h ≤ w and lattice points are ordered from top to bottom in each column and columns from left to right. The complexity in time of the algorithm is then exponential in the number of rows and linear in the number of column. The algorithm of [START_REF] Reeves | Efficient recursions for general factorisable models[END_REF] was extended in [START_REF] Friel | Recursive computing and simulation-free inference for general factorizable models[END_REF] to also allow exact draws from π(x | ψ, G ) for small enough lattices. The reader can find below an example of implementation for the general Potts model.

Following the work of [START_REF] Friel | Recursive computing and simulation-free inference for general factorizable models[END_REF], the R-package GiRaF available on CRAN proposes amongst other tools an ingenious implementation of those recursions for the autologistic, the Ising and the Potts model. Indeed, a naive implementation of the aforementioned recursions for such models can substantially increase the cost in time and memory allocation. For instance, the unnormalized Potts distribution associated to energy function (2.4) writes as

q (x | ψ, G 4 ) = n-h i=1 q i (x i:i+h | ψ, G 4 ) ,
where • for all lattice point i except the ones on the last row or last column

q i (x i:i+h | ψ, G 4 ) = exp K-1 k=0 α k 1{x i = k} + β 0 1{x i = x i+1 } + β 1 1{x i = x i+h } . (4.1)
• When lattice point i is on the last row x i+1 drops out ot (4.1), that is

q i (x i:i+h | ψ, G 4 ) = exp K-1 k=0 α k 1{x i = k} + β 1 1{x i = x i+h } . (4.2)
• The last factor takes into account all potentials within the last column

q n-h (x n-h:n | ψ, G 4 ) = exp n i=n-h K-1 k=0 α k 1{x i = k} + β 1 1{x n-h = x n } + β 0 n i=n-h+1 1{x i = x i+1 } .
One shall remark that for a homogeneous random field (i.e., parameters are independent of the location of the sites), factors (4.1) and ( 4.2) only depend on the value of the random variables X i:i+h but not on the actual position of the sites. Hence the number of factors to be computed is 2K h instead of h(w -1)K h . In term of implementation that also means factors can be computed for the different possible configurations once upstream the recursion. Furthermore with a first order neighbourhood, factor at a site merely involves its neighbour below and on its right, thereby reducing the number of possible factor to K 3 + K 2 . Finally, mention it is straightforward to extend this algorithm to hidden Markov random field since as already mention in Section 2.4 the noise corresponds to a non homogeneous potential on singleton and hence the model still writes as a general factorisable model.

Pseudo-model and variational approaches

A point of view to overcome model's bottleneck is to replace the true model with another pseudomodel selected among a collection of much simpler probability distribution, like for variational Bayes [START_REF] Jaakkola | Bayesian parameter estimation via variational methods[END_REF], or with an easily-normalised full conditional distributions [START_REF] Lindsay | Composite likelihood methods[END_REF]. Both options have been explored in the literature and we discuss some of the solutions below.

Composite likelihood

A composite likelihood [START_REF] Lindsay | Composite likelihood methods[END_REF] approximates the joint distribution as the product of easy normalised full-conditional distributions

f CL (x | θ, G ) = C i=1 f x A(i) x B(i) , θ, G wi , (5.1)
where A(i) C i=1 ⊆ P(S ) and B(i) C i=1 ⊆ P(S ) denote sets of subset of S . It has encountered considerable interests in the statistics literature and the reader may refer to [START_REF] Varin | An overview of composite likelihood methods[END_REF] for a comprehensive overview. One of the earliest approach using composite likelihood is the pseudolikelihood [START_REF] Besag | Statistical Analysis of Non-Lattice Data[END_REF] which approximates the joint distribution of x as the product of full-conditional distributions for each site i,

f pseudo (x | θ, G ) = n i=1 f (x i | x -i , θ, G ) .
(5.2)

The pseudolikelihood (5.2) is not a genuine probability distribution, except if the random variables X i are independent. The Markov property ensures that each term in the product only involves nearest neighbours, and so the normalising constant of each full-conditional is straightforward to compute. It is straightforward to prove that a unique maximum exists and it is easy to compute. [START_REF] Geman | Markov Random Field Image Models and Their Applications to Computer Vision[END_REF] demonstrate the consistency of the maximum pseudolikelihood estimator θMPLE = arg max

θ log f pseudo (x | θ, G ) .
when the lattice size tends to infinity for discrete Markov random field but maximum pseudolikelihood estimator has generally larger asymptotic variance than maximum likelihood estimator and does not achieve the Cramer-Rao lower bound [START_REF] Lindsay | Composite likelihood methods[END_REF]. Practically, this approximation has been shown to lead to unreliable estimates of θ (e.g., [START_REF] Rydén | Computational Bayesian analysis of hidden Markov models[END_REF][START_REF] Friel | Likelihood Estimation and Inference for the Autologistic Model[END_REF][START_REF] Cucala | A Bayesian Reassessment of Nearest-Neighbor Classification[END_REF][START_REF] Friel | Bayesian Inference in Hidden Markov Random Fields for Binary Data Defined on Large Lattices[END_REF].

Variational approaches and parameter estimation

Variational methods refer to a class of deterministic approaches. They consist in introducing a variational function as an approximation to the likelihood in order to solve a simplified version of an optimization problem. It has long-standing antecedents in statistical mechanics when one aims at predicting the response to the system to a change in the Hamiltonian. One important technique is based on a variational approach as the minimizer of the free energy, sometimes referred to as variational or Gibbs free energy and defined with the Kullback-Leibler divergence between a probability distribution P and the target distribution f (• | θ, G ) as

F (P) = -log Z (θ, G ) + KL (P f (• | θ, G )) .
(5.3)

Although the Kullback-Leibler divergence is not a true metric, it has the non-negative property with divergence zero if and only if distributions are equal almost everywhere. The free energy has then an optimal lower bound achieved for

P = f (• | θ, G ).
Minimizing the free energy with respect to the set of probability distribution on X allows to recover the Gibbs distribution but presents the same computational intractability. A solution is to minimize the Kullback-Leibler divergence over a restricted class of tractable probability distribution on X . This is the basis of mean field approaches that aim at minimizing the Kullback-Leibler divergence over the set of probability functions that factorize on sites of the lattice. namely for all x in X = i∈S X i ,

P(x) = i∈S P i (x i ), where P i ∈ M + 1 (X i ) and P ∈ M + 1 (X ).
The minimization of (5.3) over this set leads to fixed point equations for each marginal of P and the optimal solution is the so called mean field approximation (see for example [START_REF] Jordan | An Introduction to Variational Methods for Graphical Models[END_REF].

Variational EM algorithm

In what follows we focus on the issue of estimating the parameters of a hidden Markov random fields. Let consider y a noisy version of a Markov random field x. One can write the log-likelihood as follows

log f (y | ψ) = x∈X log f (x, y | ψ, G ) P(x) P(x)
=F (P,ψ)

+ x∈X log P(x) f (x | y, ψ, G ) P(x) =KL(P f (•|y,ψ,G ))
.

(5.4)

Solutions based on variational approaches focus on minimising the KL-term or equivalently maximising the F . For instance, this relaxation of the original issue has shown good performances for approximating the maximum likelihood estimate [START_REF] Celeux | EM procedures using mean field-like approximations for Markov model-based image segmentation[END_REF], as well as for Bayesian inference on hidden Markov random fields [START_REF] Mcgrory | Variational Bayes for estimating the parameters of a hidden Potts model[END_REF]. We discuss it below. [START_REF] Celeux | EM procedures using mean field-like approximations for Markov model-based image segmentation[END_REF] explore the opportunity to use the Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF]. To address the problem of finding the maximum likelihood estimator for the statistical model f (y | ψ), the EM algorithm iterates between two steps:

1. E step: one computes the conditional expectation of the log-joint distribution with respect to the distribution of X given Y = y at the current parameter value ψ (t)

Q(ψ | ψ (t) ) = E ψ (t) {log f (X, y | ψ, G ) | Y = y} = x∈X f (x | y, ψ (t) , G ) log f (x, y | ψ, G ) .
2. M step: one maximises Q with respect to ψ, (t) .

ψ (t+1) = arg max ψ Q ψ ψ
We refer the reader to [START_REF] Wu | On the convergence properties of the EM algorithm[END_REF] for convergence results. The EM scheme cannot be applied directly to hidden Markov random fields as it yields analytically intractable updates. The function Q can be written as

Q(ψ | ψ (t) ) = E ψ (t) {log g φ (y | X) | Y = y, } + E ψ (t) {log f (X | θ, G ) | Y = y} .
and hence requires evaluating f (• | θ, G ) and f (• | y, ψ, G ) (for computing E ψ ) which are both unavailable. Many stochastic or deterministic schemes have been proposed to handle intractability in the EM steps, such as the Gibbsian-EM [START_REF] Chalmond | An iterative Gibbsian technique for reconstruction of m-ary images[END_REF], the Monte-Carlo EM [START_REF] Wei | A Monte Carlo Implementation of the EM Algorithm and the Poor Man's Data Augmentation Algorithms[END_REF] or the Restoration-Maximization algorithm [START_REF] Qian | Estimation of parameters in hidden Markov models[END_REF]. The aim of the variational EM (VEM) is to maximize the function F instead of Q in order to get a tractable version of the EM algorithm. This shift in the formulation leads to an alternating optimization procedure which can be described as follows: let denote D a set of probability distributions on the latent space , given a current value (P (t) , θ (t) ) in D × Θ, updates with (t) , G ) , (5.5)

P (t+1) = arg max P∈D F (P, ψ (t) ) = arg min P∈D KL P f (• | y, ψ
ψ (t+1) = arg max ψ F P (t+1) , ψ = arg max θ x∈X P (t+1) (x) log π (x, y | ψ, G ) .
(5.6)

Optimising F over the class of independent probability distributions P that factorize on sites leads to mean field approximation. Generalizing an idea originally introduced by Zhang (1992), [START_REF] Celeux | EM procedures using mean field-like approximations for Markov model-based image segmentation[END_REF] have designed a class of VEM-like algorithm that uses mean field-like approximations for both f (• | y, ψ, G ) and f (• | θ, G ). To put it in simple terms mean field-like approximations refer to distributions for which neighbours of site i are set to constants. Given a configuration x in X , the Gibbs distribution f (• | θ, G ) is replaced by

P MF-like x (x | θ, G ) = i∈S P x i ; xN (i) , θ, G .
The main difference with the pseudolikelihood (5.2) is that neighbours are not random anymore and setting them to constant values leads to a system of independent variables. From this approximation, the EM path is set up with the corresponding joint distribution approximation

P MF-like (x, y | ψ, G ) = i∈S g φ (y i | x i ) P x i ; xN (i) , θ, G .
Note that this general procedure corresponds to the so-called point-pseudo-likelihood EM algorithm proposed by [START_REF] Qian | Estimation of parameters in hidden Markov models[END_REF]. The flexibility of the approach proposed by [START_REF] Celeux | EM procedures using mean field-like approximations for Markov model-based image segmentation[END_REF] lies in the choice of the configuration x that is not necessarily a valid configuration for the model. We refer the reader to [START_REF] Celeux | EM procedures using mean field-like approximations for Markov model-based image segmentation[END_REF] for further details. When the neighbours X N (i) are fixed to their mean value, or more precisely x is set to the mean field estimate of the complete conditional distribution f (x | y, ψ, G ), this results in the Mean Field algorithm of [START_REF] Zhang | The mean field theory in EM procedures for Markov random fields[END_REF]. In practice, [START_REF] Celeux | EM procedures using mean field-like approximations for Markov model-based image segmentation[END_REF] obtain better performances with their so-called Simulated Field algorithm (see Algorithm 1). In this stochastic version of the EM-like procedure, x is a realization drawn from the conditional distribution f (• | y, ψ (t) , G ) for the current value of the parameter ψ (t) . The latter is preferred to other methods when dealing with maximum-likelihood estimation for hidden Markov random field. This extension of VEM algorithms suffers from a lack of theoretical support due to the propagation of the approximation to the Gibbs distribution f (• | θG ). One might advocate in favour of the Monte-Carlo VEM algorithm of [START_REF] Forbes | Combining Monte Carlo and Mean Field-Like Methods for Inference in Hidden Markov Random Fields[END_REF] for which convergence results are available. However the Simulated Field algorithm provides better results for the estimation of the spatial parameter, as illustrated in [START_REF] Forbes | Combining Monte Carlo and Mean Field-Like Methods for Inference in Hidden Markov Random Fields[END_REF].

Approximating model choice criteria

Various approximations have been proposed to overcome intractability of (3.3) but a commonly used one, if only for its simplicity, is the Bayesian Information Criterion (BIC) that is an asymptotic estimate of the evidence based on the Laplace method for integrals (Schwarz, 1978, Kass and[START_REF] Kass | Bayes factors[END_REF]. The criterion is a simple penalized function of the maximized log-likelihood Initialization: start from an initial guess ψ (0) = θ (0) , φ (0) ; for t ← 1 to T do neighbourhood restoration: draw x(t) from π • y, ψ (t-1) , G ; E-step: compute (t) , φ (t) where

BIC(m) = -2 log f m y ψMLE + d m log(n) ≈ -
Q 1 (φ) := i∈S xi P x i ; x(t) N (i) , y i ,ψ (t-1) , G log g φ (y i | x i ) ; Q 2 (θ) := i∈S xi P x i ; x(t) N (i) , y i ,ψ (t-1) , G log P x i ; x(t) N (i) , θ, G ; M-step: set ψ (t) = θ
θ (t) = arg max θ Q 1 (θ) and ψ (t) = arg max ψ Q 2 (ψ);
end return ψ (T ) = θ (T ) , φ (T ) increases with the complexity of the model and the model with the highest posterior probability is the one that minimizes BIC. The criterion is closely related to the Akaike Information Criterion (AIC, [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF] that solely differs in the penalization term. We refer the reader to [START_REF] Kass | Bayes factors[END_REF] and the references therein for a more detailed discussion on AIC and for instance to [START_REF] Burnham | Model selection and multimodel inference: a practical information-theoretic approach[END_REF] for comparison between AIC and BIC.

BIC can be defined beside the special case of independent random variables. In the latter case the number of free parameter is, in general, not equal to the dimension of the parameter space as for the independent case. The consistency of BIC has been proven in various situations such as independent and identically distributed processes from the exponential families [START_REF] Haughton | On the Choice of a Model to Fit Data from an Exponential Family[END_REF], mixture models [START_REF] Keribin | Consistent Estimation of the Order of Mixture Models[END_REF] or Markov chains [START_REF] Csiszár | The consistency of the BIC Markov order estimator[END_REF][START_REF] Gassiat | Likelihood ratio inequalities with applications to various mixtures[END_REF]. When dealing with observed Markov random fields, aside from the problem of intractable likelihoods the number of free parameters in the penalty term has no simple formula. In the context of selecting a neighbourhood system, [START_REF] Csiszár | Consistent Estimation of the Basic Neighborhood of Markov Random Fields[END_REF] proposed to replace the likelihood by the pseudolikelihood (5.2) and modify the penalty term as the number of all possible configurations for the neighbouring sites. The resulting criterion is shown to be consistent as regards this model choice. Up to our knowledge such a result has not been yet derived for hidden Markov random fields.

As already mention, in the context of Markov random fields, difficulties are of two kind. Neither the maximum likelihood estimate θMLE nor the likelihood f m are available. Recall that in the hidden case f m requires to integrate a Gibbs distribution over the latent space configuration.

As regards the simplest case of observed Markov random field solutions have been brought by penalized pseudolikelihood [START_REF] Ji | A consistent model selection procedure for Markov random fields based on penalized pseudolikelihood[END_REF] or MCMC approximation of BIC [START_REF] Seymour | Approximate Bayes model selection procedures for Gibbs-Markov random fields[END_REF]. Over the past decade, only few works have addressed the model choice issue for hidden Markov random field from that BIC perspective. Below, we describe solutions based on pseudo-models but other attempts based on simulations techniques have been investigated [START_REF] Newton | Approximate Bayesian inference with the weighted likelihood bootstrap[END_REF].

The central question is the evaluation of the integrated likelihood (3.5). A convenient way to circumvent the issues of computing BIC is to replace the Gibbs distribution by tractable surrogates. Consider a partition of S into subsets of neighbouring sites, namely

S = C =1 A( ),
and denote by D the class of independent probability distributions P that factorize with respect to this partition, that is if X A( ) stands for the configuration space of the subset A( ), for all x in X

P(x) = C =1 P x A( ) , where P ∈ M + 1 X A( ) and P ∈ M + 1 (X ).
The most straightforward approach is to look for mean-field like approximations which take the following form

P x(x | θ, G ) = n i=1 P(x i ; xN (i) , θ, G ) ≈ f (x | θ, G ).
Put in other words, the latter is a surrogate in the class of independent probability distributions that factorize with respect to the nodes of G and for which neighbourhood of each site has been set to a constant field x. One shall remark we recover the mean-field approximations presented in Section 5.2 when x is set to the mean-field realisation. The integrated likelihood corresponding to P x is of the form

P MF-like m (y | ψ m = (θ m , φ m )) = i∈S xi g φ m (y i | x i ) P(x i ; xN (i) , θ m , G ).
This results in the following approximation of BIC BIC MF-like (m) = -2 log P MF-like m y ψMLE + d m log(n).

(5.8)

This approach includes the Pseudolikelihood Information Criterion (PLIC) of [START_REF] Stanford | Approximate Bayes factors for image segmentation: The pseudolikelihood information criterion (PLIC). Pattern Analysis and Machine Intelligence[END_REF] as well as the mean field-like approximations of BIC proposed by [START_REF] Forbes | Hidden Markov random field model selection criteria based on mean field-like approximations[END_REF]. The main difference between these criterion lies in the estimation ψ MLE and the choice of x. The idea of [START_REF] Stanford | Approximate Bayes factors for image segmentation: The pseudolikelihood information criterion (PLIC). Pattern Analysis and Machine Intelligence[END_REF] is to consider as x a configuration close to the Iterated Conditional Modes (ICM, [START_REF] Besag | On the Statistical Analysis of Dirty Pictures[END_REF]) estimate of x. In its unsupervised version ICM alternates between a restoration step of the latent states and an estimation step of the parameter ψ. PLIC is the approximation of BIC based on the output of ICM algorithm x ICM , ψ ICM . The solution proposed by [START_REF] Forbes | Hidden Markov random field model selection criteria based on mean field-like approximations[END_REF] is to use for (x, ψMLE ) the output of the VEM-like algorithm based on the mean-field like approximations of [START_REF] Celeux | EM procedures using mean field-like approximations for Markov model-based image segmentation[END_REF]. More precisely, as regards neighbourhood restoration step, they advocate in favor of the simulated field algorithm (see Algorithm 1). [START_REF] Stoehr | Hidden Gibbs random fields model selection using Block Likelihood Information Criterion[END_REF] have extended previous approaches to tractable approximations that factorize over larger sets of nodes, namely blocks of a rectangular lattice, by taking advantage of the general recursion implemented in the R-package GiRaF (see Section 4). They suggest to chose surrogates of the form

P (x ; x, A(1), . . . , A(C), θ) = C =1 P x A( ) ; X B( ) = xB( ) , θ, G ,
(5.9)

where B( ) is the border of A( ), see Figure 3, or the empty set. This leads to their so called Block Likelihood Information Criterion (BLIC) which approximates BIC as a summation of tractable normalising constant, namely BLIC is to a certain extent related to the RDA approximations of partition functions proposed by [START_REF] Friel | Bayesian Inference in Hidden Markov Random Fields for Binary Data Defined on Large Lattices[END_REF]. Indeed, let Z(θ, G ) and Z(ψ, G ) denote the respective normalizing constants of the latent and the conditional fields. Starting from the Bayes formula, BIC expression turns into

BLIC (m) = -2 C =1 log Z ψMLE , G , y A( ) , xB( ) log Z θMLE , G , xB( ) + d m log(|S |), (
BIC(m) = -2 log Z (ψ, G ) log Z (θ, G ) + d m log(n).
Similarly to [START_REF] Friel | Bayesian Inference in Hidden Markov Random Fields for Binary Data Defined on Large Lattices[END_REF], [START_REF] Stoehr | Hidden Gibbs random fields model selection using Block Likelihood Information Criterion[END_REF] approximate the intractable normalising constants by a product of tractable normalising constant defined on contiguous sub-lattices. Looking at the issue of estimating the partition function instead of estimating the Gibbs distribution has also been explored by [START_REF] Forbes | Hidden Markov random field model selection criteria based on mean field-like approximations[END_REF]. They propose to use a first order approximation of the partition function arising from mean field theory. [START_REF] Forbes | Hidden Markov random field model selection criteria based on mean field-like approximations[END_REF] argue that the latter is more satisfactory than BIC MF-like (m) in the sense it is based on a optimal lower bound for the normalising constants contrary to the mean field-like approximations. However that does not ensure better results as regards model selection.

Regarding the question of inferring the number of latent states, one might avocate in favor of the Integrated Completed Likelihood (ICL, [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. This opportunity has been explored by [START_REF] Cucala | Bayesian Inference on a Mixture Model With Spatial Dependence[END_REF] but their complex algorithm cannot be extended easily to other model selection problem such as choosing the dependency structure.

6 Sampling methods for Markov random fields

Sampling from Gibbs distribution

Sampling from a Gibbs distribution can be a daunting task due to the correlation structure on a high dimensional space, and standard Monte Carlo methods are impracticable except for very specific cases. In the Bayesian paradigm, Markov chain Monte Carlo (MCMC) methods have played a dominant role in dealing with such problems, the idea being to generate a Markov chain whose stationary distribution is the distribution of interest. The Ising model is one of these special cases where one can be drawn exactly from the model using coupling from the past [START_REF] Propp | Exact Sampling with Coupled Markov chains and Applications to Statistical Mechanics[END_REF]Wilson, 1996, Mira et al., 2001). Nevertheless, such perfect sampling scheme is often prohibitively expensive or impossible to carry out for other models. We describe below two popular solutions even though they introduce a bias.

Gibbs sampler

The Gibbs sampler is a highly popular MCMC algorithm in Bayesian analysis starting with the influential development of [START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF]. It can be seen as a componentwise Metropolis-Hastings algorithm [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF][START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF] where variables are updated one at a time and for which proposal distributions are the full conditionals themselves, see Algorithm 2. It is particularly well suited to Markov random field since by nature the intractable joint distribution is fully determined by the easy to compute conditional distributions.

Under the irreducibility assumption, the chain converges to the target distribution f (• | θ, G ), see for example (Geman and Geman, 1984, Theorem A). Note the order in which the components are updated in Algorithm 2 does not make much difference as long as every site is visited. Hence it can be deterministically or randomly modified, especially to avoid possible bottlenecks when visiting the configuration space. A synchronous version is nonetheless unavailable since updating the sites merely at the end of cycle t would lead to incorrect limiting distribution.

We should mention here that Gibbs sampler faces some well known difficulties when it is applied to the Ising or Potts model. The Markov chain mixes slowly, namely long range interactions require many iterations to be taken into account, such that switching the color of a large homogeneous area is of low probability even if the distribution of the colors is exchangeable. This ;

for t ← 1 to T do for i ← 1 to n do draw x (t) i from the full conditional P X (t) i x (t-1)
N (i) , θ ; end end return the configuration x (T ) peculiarity is even worse when the parameter β is above the critical value of the phase transition, the Gibbs distribution being severely multi-modal (each mode corresponding to a single color configuration). [START_REF] Liu | Peskun's theorem and a modified discrete-state gibbs sampler[END_REF] proposed a modification of the Gibbs sampler that overcome these drawbacks with a faster rate of convergence. Note also that in the context of Gaussian Markov random field some efficient algorithm have been proposed like the fast sampling procedure of [START_REF] Rue | Fast Sampling of Gaussian Markov Random Fields[END_REF].

Auxiliary variables and Swendsen-Wang algorithm

An appealing alternative to bypass slow mixing issues of the Gibbs sampler is the Swendsen-Wang algorithm [START_REF] Swendsen | Nonuniversal critical dynamics in Monte Carlo simulations[END_REF] originally designed to speed up simulation of Potts model close to the phase transition. Swendsen-Wang algorithm iterates two steps : a clustering step and a swapping step (see Algorithm 3), in order to incorporate simultaneous updates of large homogeneous regions (e.g., [START_REF] Besag | Spatial Statistics and Bayesian Computation[END_REF]. The clustering step relies on auxiliary random variables which aim at decoupling the complex dependence structure between the component of x and yield a partition of sites into single-valued clusters or connected components. The method set binary (0-1) conditionally independent auxiliary variables U ij which satisfy

P (U ij = 1 | x) = 1 -exp (β ij 1{x i = x j }) = p ij if i G ∼ j, 0 otherwise
with β ij ≥ 0 so that p ij takes value between 0 and 1. The latter then represents the probability to keep an egde between neighbouring sites in G . Auxiliary variables U ij induce on x a subgraph Γ(G , x) of the dependency graph G , namely the undirected graph made of edges of G for which U ij = 1, see Figure 4. During the swapping step, each cluster C of Γ(G , x) is assigned to a new state k with probability P (X C = k) ∝ exp i∈C α k . For the special but important case where α = 0, new possible states are equally likely. Also for large values of β, the algorithm manages to switch colors of wide areas, achieving a better cover of the configuration space.

For the original proof of convergence, we refer the reader to [START_REF] Swendsen | Nonuniversal critical dynamics in Monte Carlo simulations[END_REF] and for further discussion see for example [START_REF] Besag | Spatial Statistics and Bayesian Computation[END_REF]. Whilst the ability to change large set of variables in one step seems to be a significant advantage, this can be marred by a slow mixing time, namely exponential in n [START_REF] Gore | The Swendsen-Wang process does not always mix rapidly[END_REF]. The mixing time of the algorithm is polynomial in n for Ising or Potts models with respect to the graphs G 4 and G 8 but only for (t) with probability P X (t) C = k ∝ exp i∈C α k ; end return the configuration x (T ) small enough value of β [START_REF] Cooper | Mixing properties of the Swendsen-Wang process on classes of graphs[END_REF]. This was proved independently by [START_REF] Huber | A bounding chain for Swendsen-Wang[END_REF] who also derive a diagnostic tool for the convergence of the algorithm to its invariant distribution, namely using a coupling from the past procedure.

The algorithm can be extended to other Markov random field or models (e.g., [START_REF] Edwards | Generalization of the Fortuin-Kasteleyn-Swendsen-Wang representation and Monte Carlo algorithm[END_REF][START_REF] Wolff | Collective Monte Carlo updating for spin systems[END_REF][START_REF] Higdon | Auxiliary variable methods for Markov chain Monte Carlo with applications[END_REF][START_REF] Barbu | Generalizing Swendsen-Wang to sampling arbitrary posterior probabilities[END_REF] but is then not necessarily efficient. In particular, it is not well suited for latent process. The bound to the data corresponds to a non-homogeneous external field that slows down the computation since the clustering step does not make a use of the data. A solution that might be effective is the partial decoupling of [START_REF] Higdon | Discussion on the Meeting on the Gibbs Sampler and Other Markov Chain Monte Carlo Methods[END_REF][START_REF] Higdon | Auxiliary variable methods for Markov chain Monte Carlo with applications[END_REF]. More recently, [START_REF] Barbu | Generalizing Swendsen-Wang to sampling arbitrary posterior probabilities[END_REF] make a move from the data augmentation interpretation to a Metropolis-Hastings perspective in order to generalize the algorithm to arbitrary probabilities on graphs. Nevertheless, it is not straightforward to bound the Markov chain of such modifications and mixing properties are still an open question despite good results in numerical experiments.

Monte Carlo maximum likelihood estimator

The use of Monte-Carlo techniques in preference to pseudolikelihood to compute maximum likelihood estimates has been especially highlighted by [START_REF] Geyer | Constrained Monte Carlo Maximum Likelihood for Dependent Data[END_REF]. Let assume V θ,G is continuously differentiable on Θ. The gradient of the log-density (2.2) can be written as

∇ θ log f (x | θ, G ) = -∇ θ V θ,G (x) -∇ θ log Z(θ, G ).
Despite closed form is out of reach, forward-simulations from the likelihood taken at each leapfrog step can be used to provide a Monte Carlo estimate of the gradient, using the following identity,

∇ θ log Z(θ, G ) = 1 Z(θ, G ) ∇ θ Z(θ, G ) = 1 Z(θ, G ) ∇ θ X exp {V θ,G (x)} µ(dx) = X ∇ θ V θ,G (x) exp {V θ,G (x)} Z(θ, G ) µ(dx) = E θ {∇ θ V θ,G (X)} . (6.1)
So far, we have only assumed that V θ,G is continuously differentiable on Θ. However this identity holds under regularity conditions which allow one to switch the derivative and integral operators (the domain X of X is assumed to be independent of θ) and under the assumption that ∇ θ V θ,G (X) is integrable with respect to f (x | θ, G )µ(dx). When the function V θ linearly depends on the vector of parameters θ, that is

V θ,G (x) = -θ T S(x),
where S(x) = (s 1 (x), . . . , s d (x)) is a vector of sufficient statistics, expectation (6.1) is simply the first moment of the statistics with respect to f (x | θ, G )µ(dx). It is also possible to show that the log-density is concave as its Hessian matrix is the second moment of the statistics with respect to f (x | θ, G )µ(dx):

∇ 2 θ log f (x | θ, G ) = -Cov θ {S(X)} .
The maximum likelihood estimator θMLE is then the unique zero of the score function ∇ θ log f (x | θ, G ) and it satisfies

S(x) -E θMLE {S(X)} = 0. (6.2)
Hence a solution to solve problem (3.1) is to resort to stochastic approximations on the basis of equation (6.2) (e.g., [START_REF] Younes | Estimation and annealing for Gibbsian fields[END_REF][START_REF] Descombes | Estimation of Markov random field prior parameters using Markov chain Monte Carlo maximum likelihood[END_REF]. [START_REF] Younes | Estimation and annealing for Gibbsian fields[END_REF] sets a stochastic gradient algorithm converging under mild conditions. At each iteration t, the algorithm takes the direction of the gradient estimated by the value of the statistic function for one realisation of the random field x (t) , namely

θ (t+1) = θ (t) + δ t + 1 S(x) -S x (t) ,
where δ is a user-defined threshold. [START_REF] Younes | Estimation and annealing for Gibbsian fields[END_REF] provides theoretical conditions on δ which guarantee the convergence of the algorithm. However such theoretical value for δ yields in practice a step size too small to ensure the convergence to be achieved in reasonable amount of time. This can be overcome by controlling the probability of non-convergence, see the original paper by [START_REF] Younes | Estimation and annealing for Gibbsian fields[END_REF] for discussion. Another approach to compute the maximum likelihood estimation is to use direct Monte Carlo calculation of the likelihood such as the MCMC algorithm of [START_REF] Geyer | Constrained Monte Carlo Maximum Likelihood for Dependent Data[END_REF]. The convergence in probability of the latter toward the maximum likelihood estimator is proven for a wide range of models including Markov random fields. Following that work, [START_REF] Descombes | Estimation of Markov random field prior parameters using Markov chain Monte Carlo maximum likelihood[END_REF] derive also a stochastic algorithm that, as opposed to [START_REF] Younes | Estimation and annealing for Gibbsian fields[END_REF], takes into account the distance to the maximum likelihood estimator using importance sampling.

Computing posterior distributions

Markov Chain Monte Carlo allow to asymptotically sample from analytically intractable posterior distribution π. It provides a very general framework to allow estimation of functionals of the form

Θ g(θ)π(dθ),
for some function g by generating a Markov chain (θ n ) n∈N with transition kernel P which leaves π invariant. The empirical distribution so obtained leads to the following approximation

Θ g(θ)π(dθ) ≈ 1 N N n=1 g (θ i ) .
While some Bayesian estimators can be efficiently estimated with such methods via the empirical distribution, the intractability of the likelihood model in (2.2) implies, in particular, that the standard MCMC toolbox is infeasible. Indeed, we are concerned with the situation where the unnormalised posterior distribution, the right-hand-side of (3.2) is intractable. This complication results in what is often termed a doubly-intractable posterior distribution, since the posterior distribution itself is normalised by the evidence (or marginal likelihood) which is typically also intractable. For instance, proposing to move from (θ) to (θ ) in a standard Metropolis-Hastings requires the computation of the unknown normalising constants, Z(θ, G ) and Z(θ , G ),

ρ(θ, θ ) = 1 ∧ Z (θ, G ) Z θ , G q x θ , G ν(θ | θ )π(θ ) q (x | θ, G ) ν(θ | θ)π(θ) ,
where ν(θ | θ) denotes the proposal distribution to move from θ to θ . A solution, while being time consuming, is to estimate the ratio of the partition functions using path sampling [START_REF] Gelman | Simulating normalizing constants: From importance sampling to bridge sampling to path sampling[END_REF]. Starting from equation (6.1), the path sampling identity writes as log

Z (θ 0 , G ) Z (θ 1 , G ) = θ1 θ0 E θ {∇ θ V θ,G (X)}dθ.
Hence the ratio of the two normalizing constants can be evaluated with numerical integration. For practical purpose, this approach can barely be recommended within a Metropolis-Hastings scheme since each iteration would require to compute a new ratio.

The exchange algorithm [START_REF] Murray | MCMC for doubly-intractable distributions[END_REF] is a popular MCMC methods to allow sampling from doubly-intractable distributions. The exchange algorithm samples from an augmented distribution

π(θ , θ, u | x) ∝ f (x | θ, G )π(θ)ν(θ | θ)f (u | θ , G ).
whose marginal distribution in θ is the posterior distribution of interest. It extends an idea introduced by [START_REF] Møller | An efficient Markov chain Monte Carlo method for distributions with intractable normalising constants[END_REF]. The proposal of [START_REF] Møller | An efficient Markov chain Monte Carlo method for distributions with intractable normalising constants[END_REF] consists in including an auxiliary variable U whose density is the intractable likelihood itself. It follows a method based on single point importance sampling approximations of the partition functions Z(θ, G ) and Z(θ , G ). [START_REF] Murray | MCMC for doubly-intractable distributions[END_REF] develop this work further by directly estimating the ratio Z(θ, G )/Z(θ , G ) instead of using previous single point estimates. This leads to a clever algorithm to sample from the above augmented distribution, where it turns out that the ratio of intractable normalising constants drops out of the acceptance probability [START_REF] Murray | MCMC for doubly-intractable distributions[END_REF] point out that the fraction q(u | θ, G )/q(u | θ , G ) which appears above, can be considered as an single sample importance estimator of Z(θ, G )/Z(θ , G ) since it holds that

ρ θ, θ , u = 1 ∧ Z(θ, G ) Z(θ , G ) Z(θ , G )q(u | θ, G ) q(u | θ , G ) Z(θ, G ) q(x | θ , G )ν(θ | θ )π(θ ) q(x | θ, G )ν(θ | θ)π(θ)
E θ q(U | θ, G ) q(U | θ , G ) = Z(θ) Z(θ ) , (6.3)
where E θ is the expectation with respect to U ∼ f (• | θ , G ). In fact [START_REF] Alquier | Noisy Monte Carlo: convergence of Markov chains with approximate transition kernels[END_REF], consider a generalised exchange algorithm based, at each step of the algorithm, on an improved unbiased estimate of Z(θ)/Z(θ ) including multiple auxiliary draws with respect to the proposed parameter, namely,

Z (θ, G ) Z θ , G = 1 N N n=1 q(u (n) | θ, G ) q(u (n) | θ , G ) , (6.4)
where the auxiliary variables u (1) , . . . , u (N ) are drawn from f (• | θ , G ). However this socalled noisy exchange algorithm no longer leaves the target distribution invariant, nevertheless it is possible to provide convergence guarantees that the resulting Markov chain is close in some sense to the target distribution. Following the argument by [START_REF] Everitt | Bayesian Parameter Estimation for Latent Markov Random Fields and Social Networks[END_REF], any algorithm producing an unbiased estimate of the normalizing constant can thus be used in place of the importance sampling approximation and will lead to a valid procedure. An alternative to previous methods presented but neglected so far in the literature is Russian Roulette sampling [START_REF] Lyne | On Russian Roulette Estimates for Bayesian Inference with Doubly-Intractable Likelihoods[END_REF] which can be used to get an unbiased estimate of 1/Z(θ). The idea to apply MCMC

Algorithm 4: Exchange algorithm

Input: an initial guess (θ (0) , θ (0) , u (0) ), a number of iterations T Output: a sample drawn from the augmented distribution π θ, θ , u x for t ← 1 to T do draw θ from ν(• | θ (t-1) );

draw u from f (• | θ, G ); compute the Metropolis-Hastings acceptance ratio ρ θ (t-1) , θ , u = 1 ∧ q(u | θ (t-1) , G ) q(u | θ , G ) q(x | θ , G )ν(θ (t-1) | θ )π(θ ) q(x | θ (t-1) , G )ν(θ | θ (t-1) )π(θ)
Exchange move: set (θ (t) , θ (t) , u (t) ) = (θ , θ (t-1) , u) with probability ρ(θ (t-1) , θ , u);

Otherwise set (θ (t) , θ (t) , u (t) ) = (θ (t-1) , θ (t-1) , u (t-1) ); end return θ (t) , θ (t) , u (t) T t=1 methods to situation where the target distribution can be estimated without bias by using an auxiliary variable construction has appeared in the generalized importance Metropolis-Hasting of Beaumont (2003) and has then been extented by [START_REF] Andrieu | The Pseudo-Marginal Approach for Efficient Monte Carlo Computations[END_REF]. This brings another justification to the aforementioned methods and possible improvement from sequential Monte Carlo literature (e.g., [START_REF] Andrieu | Particle Markov Chain Monte Carlo methods[END_REF].

While the motivation is quite similar, the exchange algorithm is more convenient to implement whilst outperforming the single auxiliary variable method (SAVM) of [START_REF] Møller | An efficient Markov chain Monte Carlo method for distributions with intractable normalising constants[END_REF]. Indeed, SAVM requires to choose the conditional distribution q(u | x, θ) which makes it difficult to calibrate. For instance, [START_REF] Cucala | A Bayesian Reassessment of Nearest-Neighbor Classification[END_REF] stress out that a suitable choice is paramount and may significantly affect the performances of the algorithm.

A practical difficulty of implementing the exchange algorithm is the requirement to sample u from Gibbs distribution f (• | θ, G ) to guarantee a valid MCMC scheme. In all generality, it is impossible or prohibitively expensive to carry out such perfect sampling. [START_REF] Everitt | Bayesian Parameter Estimation for Latent Markov Random Fields and Social Networks[END_REF] has provided convergence results when one uses instead the final draw from a Gibbs sampler with stationary distribution f (• | θ, G ) as an approximate realisation. [START_REF] Everitt | Bayesian Parameter Estimation for Latent Markov Random Fields and Social Networks[END_REF] has notably pointed out that solely few iterations of the sampler are necessary. This approach has shown good performances in practice (e.g., , Cucala et al., 2009, Caimo and[START_REF] Caimo | Bayesian inference for exponential random graph models[END_REF].

While the application of the exchange algorithm is straightforward for Bayesian parameter inference of a fully observed Markov random field, some work have been devoted to the use of the exchange algorithm for hidden Markov random fields such as the exchange marginal particle MCMC of [START_REF] Everitt | Bayesian Parameter Estimation for Latent Markov Random Fields and Social Networks[END_REF] or the estimation procedure in [START_REF] Cucala | Bayesian Inference on a Mixture Model With Spatial Dependence[END_REF]. Though these methods produce accurate results they inherit the drawback of the exchange algorithm.

ABC model selection

Approximate Bayesian computation (ABC) has generated much activity in the literature recently as it offers a way to circumvent the difficulties of models which are intractable but can be simulated from. We refer the reader to [START_REF] Marin | Approximate Bayesian Computational methods[END_REF] and the references therein for a comprehensive overview on the method. When performing parameter estimation, the method is particularly well suited for problems where the likelihood function does not admit an algebraic form, a situation where MCMC methods are at a loss. We believe that the benefit of ABC for parameter estimation of a Gibbs random field is questionable. Hence we solely focus on model selection in this part.

To approximate m MAP , ABC starts by simulating numerous triplets (m, θ m , y) from the joint Bayesian model, see Algorithm 5. Afterwards, the algorithm mimics the Bayes classifier (3.4): it approximates the posterior probabilities by the frequency of each model number associated with simulated y's in a neighbourhood of y obs . In a metric space (Y , ρ), this neighbourhood is define by the ball B( , y obs ) of radius centered at y obs . If required, we can eventually predict the best model with the most frequent model in the neighbourhood, or, in other words, take the final decision by plugging in (3.4) the approximations of the posterior probabilities. However such a naive implementation is typically infeasible as the data usually lies in a space of high dimension and the algorithm faces the curse of dimensionality. Put in other words, sample dataset in the neighbourhood of y occurs with an prohibitively low probability. The ABC algorithm performs therefore a (non linear) projection of the observed and simulated datasets onto some Euclidean space of reasonable dimension d via a function s, composed of summary statistics S(•) = s 1 (•), . . . , s M (•) that are the concatenation of the summary statistics of each models with cancellation of possible replicates.. The neighbourhood of y obs is thus defined as simulations whose distances to the observation measured in terms of summary statistics, i.e., ρ S(y), S(y obs ) , fall below a threshold . The accepted particles (m (t) , y (t) ) at the end of Algorithm 5 are distributed according to the pseudo-target π(m | ρ S(y), S(y obs ) ≤ ) and the estimate of the posterior model distribution is given by π m y obs = 1 m (t) = m, ρ S y (t) , S y obs ≤ 1 ρ S y (t) , S (y obs ) ≤ .

ABC hence presents two level of approximations arising from the size of the neighbourhood and the introduction of summary statistics.

The choice of such summary statistics presents major difficulties that have been especially highlighted for model choice [START_REF] Robert | Lack of confidence in approximate Bayesian computation model choice[END_REF][START_REF] Didelot | Likelihood-free estimation of model evidence[END_REF]. When the summary statistics are not sufficient for the model choice problem, [START_REF] Didelot | Likelihood-free estimation of model evidence[END_REF] and [START_REF] Robert | Lack of confidence in approximate Bayesian computation model choice[END_REF] found that the above probability can greatly differ from the genuine π(m | y obs ). Model selection between fully observed Markov random fields whose energy function is of the form V θ,G (y) = θ T s(y) is a surprising example for which ABC is consistent. Indeed [START_REF] Grelaud | ABC likelihood-free methods for model choice in Gibbs random fields[END_REF] have pointed out that vector of summary statistics S(•) = s 1 (•), . . . , s M (•) is sufficient for each model but also for the joint parameter across models (M , θ 1 , . . . , θ M ). This allows to sample exactly from the posterior model distribution when = 0. However the fact that the concatenated statistic inherits the sufficiency property from the sufficient statistics of each model is specific to exponential families [START_REF] Didelot | Likelihood-free estimation of model evidence[END_REF]. When dealing with model choice between hidden Markov random fields, we fall outside of the exponential families due to the bound to the data. Thus we face the major difficulty outlined by [START_REF] Robert | Lack of confidence in approximate Bayesian computation model choice[END_REF]: it is almost impossible to build a sufficient statistic of reasonable dimension, i.e., , of dimension much lower than the

  X A . Denote by S \ A = -A the complement of A in S . When modelling local interactions, the sites are lying on an undirected graph G which (a) (b)

Figure 1 :

 1 Figure 1: First and second order neighbourhood graphs G with corresponding cliques. (a) The four closest neighbours graph G4. neighbours of the vertex in black are represented by vertices in gray. (b) The eight closest neighbours graph G8. neighbours of the vertex in black are represented by vertices in gray. (c) Cliques of graph G4. (d) Cliques of graph G8.

Figure 2 :

 2 Figure 2: Phase transition for a 2-states Potts model with respect to the first order and second order 100×100 regular square lattices. (a) Average proportion of homogeneous pairs of neighbours. (b) Variance of the number of homogeneous pairs of neighbours.

  Model choice is a problem of probabilistic model comparison. Assume we are given a set M = {m : 1, . . . , M } of stochastic models with respective parameter spaces Ψ m embedded into Euclidean spaces of various dimensions. Bayesian approach to model selection considers the model itself as an unknown parameter of interest. The joint Bayesian distribution sets: a prior π on the model space, a prior density π m on each parameter space Ψ m with respect to a reference measure (often the Lebesgue measure of the Euclidean space), the likelihood f m of the data y within each model. On the extended parameter space Ψ = M m=1 {m} × ψ m , the Bayesian analysis targets posterior model probabilities, that is the marginal in M of the posterior distribution for (m, ψ 1

  2 log e(y | m), (5.7) where d m is the number of free parameters of model m (usually the dimension of Ψ m ) and n = |S | is the number of sites. The d m log(n) term corresponds to a penalty term which Algorithm 1: Simulated Field algorithm Input: an observation y, a number of iterations T Output: an estimate of the likelihood maximum estimator ψMLE

Figure 3 :

 3 Figure 3: Factorisation of the Gibbs distribution approximation over a set of contiguous block A( ) with border B( ) set to a constant field x.

  5.10)where Z θ, G , xB( ) is the partition function of block A( ) with fixed border xB( ) and Z θ, G , y A( ) , xB( ) is the partition function of the conditional random field X A( ) knowing Y A( ) = y A( ) and X B( ) = xB( ) . The latter has shown better performances for estimating the number of latent component K as well as for the selection of the underlying dependency structure G . The criterion nonetheless lack some theoretical support and is limited to regular lattices.

Algorithm 2 :

 2 Gibbs sampler Input: a parameter θ, a number of iterations T Output: a sample x from the joint distribution f (• | θ, G ) Initialization: draw an arbitrary configuration x (0) = x

Figure 4 :;

 4 Figure 4: Auxiliary variables and subgraph illustrations for the Swendsen-Wang algorithm. (a) Example of auxiliary variables Uij for a 2-states Potts model configuration on the first order square lattice. (b) Subgraph Γ(G4, x) of the first order lattice G4 induced by the auxiliary variables Uij.

  The 2-states Potts model is equivalent to the Ising model up to a constant for interaction parameter β, that is β Potts = 2β Ising .

Algorithm 5: ABC model choice algorithm

Input: an observation y obs , summary statistics S, a number of iterations T , an empirical quantile of the distance T Output: a sample from the approximated target of π • | S(y obs )

save m (t) , θ (t) , S y (t) ← {m, ψ, S (y)};

end sort the replicates according to the distance ρ S y (t) , S y obs ; keep the T first replicates; return the sample of accepted particles dimension of X .

Beyond the seldom situations where sufficient statistics exist and are explicitly known, [START_REF] Marin | Relevant statistics for Bayesian model choice[END_REF] provide conditions which ensure the consistency of ABC model choice but the latter are difficult, if not impossible, to check in practice. To answer the absence of sufficient statistics as well as aforementioned theoretical conditions, very few has been accomplished in the context of ABC model choice. One solution would be to rely on the approach of [START_REF] Prangle | Semi-automatic selection of summary statistics for ABC model choice[END_REF].

The statistics S(y) reconstructed by [START_REF] Prangle | Semi-automatic selection of summary statistics for ABC model choice[END_REF] have good theoretical properties (those are the posterior probabilities of the models in competition). Nevertheless the methods requires a pilot ABC run which is time consuming and can lead to poor approximations. Alternatively, [START_REF] Stoehr | Adaptive ABC model choice and geometric summary statistics for hidden Gibbs random fields[END_REF] have proposed to use summary statistics based on connected components of induced graphs in order to select between dependency structures of hidden Markov random fields. Beside, the specific result on Markov random fields, they derive an adaptive scheme to select the most informative set of statistics based on a local error rate. The main drawback of the method is that its use is limited to low dimensional vector of statistics. To overcome this issue, the summary statistics proposed by [START_REF] Stoehr | Adaptive ABC model choice and geometric summary statistics for hidden Gibbs random fields[END_REF] could be used among others within the ABC random forest procedure of [START_REF] Pudlo | Reliable abc model choice via random forests[END_REF].


R J Dikken
email: r.j.dikken@tudelft.nl

On the relation between emergent space-time and the Schrödinger equation

Keywords: holographic principle, entropy, quantum mechanics, information, emergent spacetime

This work shows the relation between emergent space-time and the Schrödinger equation. An additional inertial contribution from an effective potential field is incorporated into Newton's second law through an extended definition of the distance of a particle from a holographic screen and the associated change in entropy. From the same principles, the Schrödinger equation is derived, showing an intimate relation between emergent space-time and quantum mechanics.

I. INTRODUCTION

A long standing challenge in physics is the unification of the laws of quantum mechanics and gravity. A candidate theory that has been developed to this end is M-theory, being the collection of several superstring theories [START_REF] Rickles | A Brief History of String Theory: From Dual Models to M-Theory[END_REF][START_REF] Witten | String theory dynamics in various dimensions[END_REF], showing an intriguing complexity. One of the developments resulting from string theory is the holographic principle that states that a 3+1 dimensional world can be thought of through dimensionality reduction as a 2+1 dimensional world, i.e. the information of a 3+1D world is encoded on a 2+1D world [START_REF] Hooft | Dimensional Reduction in Quantum Gravity[END_REF][START_REF] Susskind | The world as a hologram[END_REF]. This explicitly results in the AdS/CFT correspondence, relating the anti-de-Sitter spaces including quantum gravity expressed in M-theory to conformal field theories describing elementary particles [START_REF] Maldacena | The Large N limit of superconformal field theories and supergravity[END_REF].

In black hole thermodynamics studies it has been shown that a clear relation exists between gravity and thermodynamics [START_REF] Bekenstein | Black holes and entropy[END_REF][START_REF] Bardeen | The Four laws of black hole mechanics[END_REF][START_REF] Unruh | Notes on black hole evaporation[END_REF][START_REF] Jacobson | Thermodynamics of spacetime: The Einstein equation of state[END_REF][START_REF] Strominger | Microscopic origin of the Bekenstein-Hawking entropy[END_REF]. Using similar relations it has been shown, with the holographic principle as guiding assumption, that space-time and Newton's laws of motion and gravity emerge naturally from entropic origins [START_REF] Verlinde | On the Origin of Gravity and the Laws of Newton[END_REF][START_REF] Pan | Influence on the entropic force by the virtual degree of freedom on the holographic screen[END_REF][START_REF] Chaichian | On gravity as an entropic force[END_REF][START_REF] Chaichian | On entropic gravity: The entropy postulate, entropy content of screens and relation to quantum mechanics[END_REF][START_REF] Diósi | Note on possible emergence time of Newtonian gravity[END_REF]. In line with this, the role of information associated with space-time is shown by the appearance of the bound on information flow, i.e. the quantum of thermal conductance [START_REF] Dikken | Information, uncertainty and holographic action[END_REF]. Also the Heisenberg uncertainty relation derives from the bound on information flow [START_REF] Dikken | Information, uncertainty and holographic action[END_REF]. This strongly points to a relation between thermodynamics, gravity/emergent space-time and quantum mechanics [START_REF] Padmanabhan | Thermodynamical aspects of gravity: new insights[END_REF].

In this work, the relations between quantum mechanics and thermodynamics and emergent space-time are explicitly shown. The role of inertia in Newton's second law is revisited in Section II from an entropic perspective. In Section III it is shown how the Schrödinger equation naturally emerges from the first principles used in Ref. [START_REF] Verlinde | On the Origin of Gravity and the Laws of Newton[END_REF] to derive Newton's second law and Newtonian gravity.

II. INERTIA IN NEWTON'S SECOND LAW

To derive Newton's second law of motion, it is assumed in Ref. [START_REF] Verlinde | On the Origin of Gravity and the Laws of Newton[END_REF] that a particle with mass m at a distance ∆x from a holographic screen of temperature T merges with the microscopic degrees of freedom of the screen, leading to a change in entropy ∆S. This situation is schematically depicted in Fig. 1. In Ref. [START_REF] Verlinde | On the Origin of Gravity and the Laws of Newton[END_REF] the distance is defined as the Compton wavelength ∆x = mc , where is the Planck constant and c is the speed of light. This leads to a ∆x-dependent change in entropy: ∆S = 2k B mc ∆x, where k B is the Boltzmann constant. An extension of this picture in this work is to additionally associate an effective potential U with the particle. This changes the definition of ∆x, since momentum is no longer only given by mc, but an additional contribution comes through the effective potential, p U = t 0 F U ∆t. Since the particle merges with the microscopic degrees of freedom of the screen, the same geodesics apply for the particle as for the screen. The geodesics of the holo- graphic screen are light-like and thus the screen is described by a null-metric, ∆s 2 = 0, which entails the relation ∆t = ∆x/c [START_REF] Dikken | Information, uncertainty and holographic action[END_REF][START_REF] Munkhammar | On Non-Equilibrium Thermodynamics of Space-Time and Quantum Gravity[END_REF]. Using this, the relation between the additional momentum and the effective potential is obtained as:

p U = - t 0 ∆U ∆x ∆t = - U/2 0 ∆U c = -U/(2c).
Since every bit carries an amount of energy 1 2 k B T [START_REF] Verlinde | On the Origin of Gravity and the Laws of Newton[END_REF], the integration is performed from zero to U/2. Hence, the new definition for ∆x becomes:

∆x = mc -U/(2c) . (1)
The implemented momentum definition is analogous to the modification of the momentum operator in a vector potential field, p → p -qA, commonly applied in the quantum mechanical description of particles in a vector potential field. Here we do not specify the characteristics of the effective potential U , which can be for instance an electro-static potential or a gravitational potential, etc. Following Ref. [START_REF] Verlinde | On the Origin of Gravity and the Laws of Newton[END_REF], we obtain a general linear relation between the change in entropy and the distance from the screen:

∆S = 2k B mc -U/(2c) ∆x. ( 2
)
In general, the entropic force in the emergent bulk spacetime is proportional to the change in entropy of the screen, F ∆x = T ∆S. Substituting the change in entropy of Eq. 2, the general expression of the entropic force leads to

F = 2k B T mc -U/(2c) . (3)
Through Unruhs law, the temperature T and the acceleration a are related,

k B T = 1 2 a c . (4)
As explained in Ref. [START_REF] Verlinde | On the Origin of Gravity and the Laws of Newton[END_REF], in this context this relation should be thought of as the temperature that is required for the acceleration and not as the temperature caused by the acceleration. The last step to obtain Newton's second law of motion, with an additional inertial contribution from the effective potential, is the substitution of Eq. 4 in Eq. 3:

F = m - U 2c 2 a. (5)
The derived expression is an extension of Newton's second law and states that inertia is not only related to mass, but also due to interaction with an effective potential field. The idea of additional inertia due to interaction with a potential field is long established, first by J.J. Thomson [START_REF] Thomson | On the electric and magnetic effects produced by the motion of electrified bodies[END_REF] for an electron in a magnetic field.

Based on this knowledge, the fundamental assumption in Eq. 2 is verified through the derivation of Eq. 5. This assumption is used in the following section to derive the Schrödinger equation.

III. ENTROPY AND THE SCHR ÖDINGER EQUATION

In this section, the Schrödinger equation is derived using the extension from Section II for the definition of ∆S. To this end, a relation between the entropy S and the wave-function ψ is required. This is obtained through the holographic relation between S and the effective action I [START_REF] Padmanabhan | Is gravity an intrinsically quantum phenomenon? From the entropy of space-time and the principle of equivalence[END_REF][START_REF] Caravelli | Holographic effective actions from black holes[END_REF]:

I = - S k B . ( 6
)
Important to note is that Ref. 6 states a holographic relation, i.e. the effective action I relates to the emergent bulk and the entropy S relates to the screen. Furthermore, the wave-function and the effective action are related by ψ ≡ exp(i I/ ) [START_REF] Kompaneyets | Theoretical physics[END_REF]. Hence, the relation between entropy of the screen and wave-function of the particle in the emergent bulk is given by S = i k B ln ψ. Note that this relation does not entail that the entropy is complex, since the wavefunction itself is complex. The change in entropy is given by

∆S = i k B ∆(ln ψ). (7)
Using the fact that ∆(ln ψ) = ∆ψ/ψ and equating Eq. 2 and Eq. 7, the following relation is obtained:

i ∆ψ ψ = 2 mc -U/(2c) ∆x. (8)
Rearranging the terms and multiplying both sides of the equation by i 2 ∆ ∆x , we obtain

- 2 2m ∆ 2 ψ ∆x 2 = i c ∆ψ ∆x -i U 2mc ∆ψ ∆x . (9)
For the moment we treat ∆ ∆x not as an operator, but simply as a variable, since it relates to the distance of the particle from the holographic screen. Rewriting the first term on the right-hand side of Eq. 9 by using ∆t = ∆x/c (see Section II), leads to

- 2 2m ∆ 2 ψ ∆x 2 = i ∆ψ ∆t -i U 2mc ∆ψ ∆x . ( 10
)
For the second term on the right-hand side of Eq. 10, an explicit relation for ∆ψ ∆x is obtained by rewriting Eq. 8, ∆ψ ∆x = -2 i mc-U/(2c) ψ. By substituting this relation in Eq. 10, rearranging the terms and assuming infinitesimal variation, the following expression is obtained:

- 2 2m d 2 ψ dx 2 + U - U 2 2mc 2 ψ = i d ψ dt . (11)
The obtained relation is a relativistic Schrödinger equation, where the relativistic effects are included through the term U 2 2mc 2 . If we make a connection with a real material and look at for instance the behavior of electrons in graphene, we observe that the electrons behave relativistically as quasi-particles with reduced mass. Hence, the term U 2 2mc 2 both includes relativistic behavior due to strong field interactions, where particles are accelerated to high velocities, and relativistic behavior through reduced mass.

In the weak interaction limit, U 2 2mc 2 → 0, Eq. 11 simplifies to the known Schrödinger equation:

- 2 2m d 2 ψ dx 2 + U ψ = i d ψ dt . (12)
IV. CONCLUSIONS

This work uncovers the relation between quantum mechanics and emergent space-time. It is found that there exists a direct relation between the Schrödinger equation and the entropic principles from which space-time and gravity are emergent in Verlinde's conceptual theory.

The results in this work can be viewed from three perspectives: 1) emergent space-time has a quantum mechanical origin, 2) quantum mechanics emerges from entropic origins and 3) there is a duality between the quantum mechanical description and the entropic description.

It remains an open question which perspective will prove to be the correct one.

V. ACKNOWLEDGEMENT I would like to thank prof. Erik Verlinde for useful comments on this work and dr. Marcel Sluiter for encouragement to present this work.

FIG. 1 .

 1 FIG. 1. Schematic representation of a particle with mass m in an effective potential field U that merges at distance ∆x with the microscopic degrees of freedom of a holographic screen with temperature T leading to a change in entropy ∆S.