{'fraction_non_alphanumeric': 0.0392624476500698, 'fraction_numerical': 0.007871878393051032, 'mean_word_length': 4.3950836820083685, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Inference and Optimal Design for Nearest-Neighbour Interaction Models

We consider problems of Bayesian inference for a spatial epidemic on a graph, where the final state of the epidemic corresponds to bond percolation, and where only the set or number of finally infected sites is observed. We develop appropriate Markov chain Monte Carlo algorithms, demonstrating their effectiveness, and we study problems of optimal experimental design. In particular, we demonstrate that for lattice-based processes an experiment on a sparsified lattice can yield more information on model parameters than one conducted on a complete lattice. We also prove some probabilistic results about the behaviour of estimators associated with large infected clusters.


I. INTRODUCTION
Many real-world phenomena can be modelled by random graphs, or more generally, by dynamically changing random graphs. Specifically, host-pathogen biological systems that may combine primary and nearestneighbour or long-range secondary infection processes can be efficiently described by spatio-temporal models based on random graphs evolving in time [1].
Although a continuous observation of an epidemic is not always possible, a spatial 'snapshot' may provide one with some, albeit highly incomplete, knowledge about the epidemic. In terms of the model this knowledge results in a random graph realised in some metric space. Moreover, under some circumstances it is not even possible to observe some or all of the edges of such a random graph-all one would know then are the vertices which correspond to the infected sites, i.e. to those sites which interacted as a result of the evolution of the process under consideration.
One particular application refers to the colonisation of susceptible sites, such as seeds or plants grown on a lattice, by virus, fungal, or bacterial pathogens with limited dispersal abilities. A typical example is the spread of infections through populations of seedlings by the fungal pathogen, Rhizoctonia solani Kühn. This economicallyimportant pathogen is wide spread with a remarkably wide host range [2]. In addition to its intrinsic economic importance, it has been extensively used as an experimental model system to test epidemiological hypotheses in replicated microcosms [3,4] and to study biological control of pathozone behaviour by an antagonistic fungus and disease dynamics [5]. Transmission of infection be-tween plants occurs by mycelial growth from an infected host, with preferential spread along soil surfaces-hence the missing information about the structure of interactions.
The spread of infections with limited dispersal abilities among plants can be viewed as a spatial SIR (susceptible → infective → removed) epidemic with nearest-neighbour secondary infections and removals, and can be related to a percolation process on a finite static network [1,6]. Percolation has also been widely used to model various phenomena of disordered media in physics, chemistry, biology, engineering and materials science [7,8], and also forest fires [9,10].
Bayesian estimation for percolation models of disease spread in plant populations in the context of the spread of Rhizoctonia solani has been presented in [1]. In [6] the spread of this soil-borne fungal plant pathogen among discrete sites of nutrient resource is studied using simple concepts of percolation theory; a distinction is made between invasive and noninvasive saprotrophic spread. The authors of [1,6] formulated statistical methods for fitting and testing percolation-based, spatio-temporal models that are generally applicable to biological or physical processes that evolve in time in spatially distributed populations. Estimation of spatial parameters from a single snapshot of an epidemic evolving on a discretised grid under the assumption that fundamental spatial statistics are near equilibrium is studied in [11].
The difficulties in performing inference for these models in the presence of observational uncertainty or incomplete observations can be overcome to an extent by employing a Bayesian approach and modern powerful computational techniques-mainly Markov chain Monte Carlo (McMC), e.g. [12]. McMC methods often offer important advantages over existing methods of analysis. Particularly, they allow a much greater degree of modelling flexibility, although the implementation of McMC methods can be problematic because of convergence and mixing difficulties which arise due to the amount and nature of missing data.
An aspect which has received little attention in the context of the described models is that of experimental design. Statisticians have investigated the question of experimental design in the Bayesian framework (see [13] for a review). The work of Müller and others [14,15] examined the ways of identifying designs that maximise the expectation of a utility function.
The goal of this paper is to study inference and optimal design problems for a model which is obtained by taking a final snapshot of nearest-neighbour disease spread spatio-temporal dynamics and relating it to the percolation process on a regular grid. We use the utility-based Bayesian framework for this purpose and discuss generic issues that arise in this context.
The paper is arranged as follows. In Section 2 we describe the percolation model considered and describe the partial sampling scenarios that form the focus of the paper. In Section 3 we describe how, by embedding the problem in a framework of Bayesian inference, it is possible to estimate likelihood functions for percolation models under the respective sampling scenarios using Markov chain Monte carlo (McMC) methods. Moreover we prove some theoretical results regarding the limiting properties of the likelihood for the case when only cluster size is observed, and elicit interesting connections with percolation thresholds. In Sections 5 and 6 we turn our attention to the question of experimental design and, in particular, the identification of experimental lattices that can maximise the information on the percolation probability yielded by an experimental process observed on it, when only the vertex set of the connected cluster is observed. We demonstrate that the optimal lattice may be one which includes regions which have been sparsified by removal of certain vertices and edges. These results have important implications for the design of experiments on, for example, fungal pathogen spread in populations of agricultural plants, where hosts are typically arranged on a lattice.

A. Model description
Disease spread as a result of (typically) short-range contact between, for example, plants can be modelled as a transmission process on an undirected graph. Nodes, or vertices, of the graph correspond to possible locations of plants, and edges of the graph link locations which are considered to be neighbours. In a classical SIR model each node, or vertex, of the graph is in one of three states: either it is occupied by a healthy, but susceptible, plant (state S ), or it is occupied by an infected and infectious plant (state I ), or finally it is empty, any plant at that location having died and thus being considered removed (state R). A plant at node i, once infected (or from time 0 if initially infected), remains in the infected state I for some random lifetime τ i after which it dies, so that node i then remains in the empty state R ever thereafter. During its infected lifetime the plant at node i sends further infections to each of its neighbouring nodes j as a Poisson process with rate λ ij (so that the probability that an infection travels from i to j in any small time interval of length h is λ ij h+ o(h) as h → 0 while the probability that two or more infections travel in the same interval is o(h) as h → 0); any infection arriving at node j changes the state of any healthy plant there to infected, and otherwise has no effect. All lifetimes and infection processes are considered independent. The initial state of the system is typically defined by one or more nodes being occupied by infected plants, the remaining nodes being occupied by healthy plants. The epidemic may die out at some finite time at which the set of infected nodes first becomes empty, or, on an infinite graph only, it is possible that it may continue forever.
Thus, for any infected node i, the event E ij that any neighbouring node j receives at least one infection from node i has probability p ij = 1 − E[exp(−λ ij τ i )] (where E denotes expectation). Note that, for any given node i, even though the infection processes are independent, the events E ij are themselves independent if and only if the random lifetime τ i is a constant. We now suppose that this is the case and that furthermore, for all ordered pairs (i, j) of neighbours, we have p ij = p for some probability p. Suppose further that it is possible to observe neither the time evolution of the epidemic nor the edges of the graph by which infections travel, but only the initially infected set of nodes and the set of nodes which are at some time infected and thus ultimately in the empty state R. It is then not difficult to see, and is indeed well known [16,17], that the epidemic may be probabilistically realised as an unoriented bond percolation process on the graph in which each edge is independently open with probability p (or closed with probability 1 − p), and in which the set of nodes which are ever infected consists of those nodes reachable along open paths (chains of open edges) from those initially infected. (Note that the ability to use an unoriented bond percolation process requires both the assumptions that the above events E ij are independent and that p ij = p ji for all i, j; in the absence of either of these assumptions one would in general need to consider an oriented process with the appropriate dependence structure [28].) In the present paper we consider the epidemic to take place on some subset P of Z d , where we allow P = Z d as a possibility. Two sites (nodes) are considered neighbours if and only if they are distance 1 apart and thus Z d induces a graph, called the d-dimensional cubic lattice L d = (Z d , E d ), with the set E d of edges connecting neighbours. Thus in the case P = Z 2 , for example, each node has 4 neighbours. Notice that the set P also induces a graph with respect to the lattice L d -we shall denote this graph Π. This may be considered as a model for nearest-neighbour interaction.
We assume furthermore that initially there is a single infected site, and that all other sites in Π are occupied by healthy individuals.

B. Incomplete observations
The probability p introduced above is considered to be unknown, possibly depending, through the experimental design, on other parameters (for example, on the spacing of the grid). We consider inference under each of the following two scenarios: (i) the observations consist of the set of sites which are ever infected, so that the routes by which infections travel are not observed; note that, in terms of the bond percolation process, this corresponds to one's knowledge of the connected component containing the initially infected site-the location of this site within the component not being relevant to inference for p (see below); (ii) all that is observed is the size of the set of sites which are ever infected.
We refer further to the former of these two scenarios as S1 and to the latter scenario as S2.
Since for any given set of initially infected sites, the distribution of the set of ever-infected sites in an SIR epidemic with constant infectious times is the same as for the corresponding unoriented bond percolation process (with the same initial set of the nodes involved), a final snapshot of such epidemic can be seen as an open cluster of the corresponding percolation process. Figure 1 depicts an open cluster obtained by simulation of percolation process on the integer lattice Z 2 for p = 0.478. This connected component containing the origin can be seen as a final (and finite) outbreak of an SIR epidemic process of the kind discussed above. The origin (or, indeed, any other vertex of the open cluster) may be considered to be the site where the initially inoculated individual has been placed. Clearly, the realised bond structure is not the only possible way which results in the site configuration from Figure 1. However, the distribution of this site configuration as an extinct SIR epidemic coincides with that of the corresponding unoriented bond percolation process.

III. PARAMETER ESTIMATION
A. Scenario S1 Let Π be a (proper or improper) subgraph of L d which contains the origin and let C be an open cluster of a percolation process on Π containing the origin. The set of nodes C represents a snapshot of an extinct outbreak of our spatial SIR epidemic evolved on Π ⊆ L 2 . Let us introduce some additional notions. Let G = (V, E) be a locally finite graph and let G ′ = (V ′ , E ′ ) be a subgraph of G (that is V ′ ⊆ V and E ′ ⊆ E). By the saturation of the graph G ′ with respect to G we understand the graphG = (Ṽ ,Ẽ) such that We denote the saturation of G ′ with respect to G by Satur G G ′ or, in cases when it is clear from the context with respect to what graph the saturation takes place, by Satur G ′ . A graph G ′ whose saturation coincides with itself is called a fully saturated graph. For example, the fully saturated graph (with respect to L 2 ) is obtained from the graph depicted in Figure 1 by connecting all pairs of neighbouring black sites (according to the 4neighbourhood relationship). Note that the operation of saturation may also be applied solely to a subset of vertices of the original graph, since it does not make use of the edges of the subgraph-operand.
To distinguish between the boundary points of a graph and their neighbours, which are not in the graph, we introduce the notions of the surface and the frontier of the graph (again, with respect to another graph). Let us denote by ∂G the surface of G in Π, G ⊆ Π, that is to say the set ∂G := {x ∈ G : ∃y ∈ Π\G s.t. x and y are neighbours in Π}, and by Γ G the frontier of G in Π, i.e. the set ∂(Π\G).
In order to identify the likelihood function we introduce the set G(C) of all connected subgraphs of Π with C as a vertex set. Note that the set G(C) is necessarily nonempty. For each G ∈ G(C) the number of edges between the vertices of the graph G and the elements of its frontier Γ G is constant-we denote it by w C . Finally, we denote the total number of edges in G by e(G).
The probability that C represents the set of everinfected sites and that the edges of G correspond to those routes along which the infection travelled is P p (G) = p e(G) (1 − p) e(Satur C)−e(G)+wC , and the likelihood function associated with the observed set C of ever-infected sites is given by Hence, under assumption of a uniform prior for p, its posterior distribution π(p | C) is a mixture of beta distributions: It is not feasible to calculate π(p | C) in the above form, since it is hard to enumerate all corresponding graphs and, thus, to calculate efficiently the coefficients r(k). We describe therefore an McMC algorithm that allows one to sample from the distribution π(p | C) under the uniform prior on p, that is, effectively, to estimate the likelihood function of p.
Our Markov chain explores the joint space (0, 1)×G(C) of possible values for p and graphs from G(C). The stationary distribution of the chain is the joint posterior distribution of p and G ∈ G(C). The description of the chain is given in Algorithm 1 (Appendix A 1). This Markov chain explores the set of all connected graphs G(C) by simply deleting or adding an edge from the current graph preserving the connectivity of the given site configuration C.
We apply Algorithm 1 to the site configuration C from Figure 1 (black dots only). Figure 2 depicts the likelihood function of the model parameter for the complete observations (i.e. nodes and edges of the cluster) and a histogram of the sample obtained from the posterior distribution π(p | C) when the prior distribution is uniform on the interval (0, 1).
Mixing properties of Markov chains described by Algorithm 1 are briefly discussed in Appendix B.

B. Scenario S2
Under this scenario only the size n of the outbreak of our SIR epidemic evolving on Π = L d is given. Let G n be the set of all possible connected graphs on n vertices including the origin. These graphs represent the outbreaks of the size n and we distinguish all isomorphic graphs which have different locations or orientations.
We denote the number of edges of Π between the vertices of the graph G ∈ G n and the vertices of its frontier Γ G by w(G).
Given the epidemic size n, the inference on p involves evaluation of the likelihood function L n (p) := P p (|C| = n) which can be represented as follows: As previously, under assumption of a uniform prior for p its posterior distribution π(p | |C| = n) is a mixture of beta distributions: This, again, represents a hard enumeration problem. However, inference on p can be made using the McMC technique. Algorithm 2 (Appendix A 1) contains a description of a Markov chain which serves the purpose of sampling from the posterior distribution π(p | |C| = n), given the prior distribution of p is uniform on [0,1]. The chain explores the joint space (0, 1) × G n of possible values for the percolation parameter p and all possible connected graphs on n nodes, by deleting a vertex from a graph and adding a vertex from its frontier.  Mixing properties of Markov chains described by Algorithm 2 are briefly discussed in Appendix B. Two video examples of dynamic graph updates illustrating Algorithms 1 and 2 can be found at the WEB address [18].

C. Convergence analytical results for S2
Percolation exhibits a phenomenon of criticality, this being central the percolation theory: as p increases, the sizes of clusters (connected components) also increase, and there is a critical value of p c at which there appears a cluster which dominates the rest of the pattern. If p < p c , then with probability one all open clusters are finite, but when p > p c there is a single infinite open cluster almost surely [19]. Note, that the critical percolation parameter p c depends on the dimension d of the lattice. Bond percolation on the square lattice seems to be most studied to date of all percolation processes. The critical probability p c in the case of a square lattice is 1 2 . What follows, however, holds for any lattice Z d , d ≥ 3.
Percolation theory tells us that if p < p c , then χ(p) < ∞. When p > p c , then χ(p) = ∞ and the function χ is not of a much interest in this case. Instead, one studies the function χ f (p) = E p [|C| : |C| < ∞]. The function χ(p) (χ f (p)) monotonically increases to infinity as p ↑ p c (p ↓ p c ). It is known [19] that there is no infinite open cluster when p = p c for percolation on the square lattice Z 2 and on the lattice Z d , where d ≥ 19; this is also believed to be true for all other values of d. How likely it is to observe an open cluster of size n when n is very large? What value of p should one suggest if one happened to observe a large epidemic of size n?
Intuitively, it is unlikely that very large but finite epidemics may be obtained for values of p significantly different from p c (if p ≪ p c , then it is unlikely that the epidemic will attain a very large size; if p ≫ p c , then it is equally unlikely that, having attained a sufficiently large size n, the epidemic will have burned out). Intuition suggests therefore that the likelihood function for p, given that the size n of the connected component containing the origin is growing, should be increasingly concentrated around p c .
Let L n (p) = P p (|C| = n) be the probability that an open cluster C from a percolation process with the edge density p is of size n. (Equivalently, assume that we observe an SIR spread of infection via nearest-neighbour interactions on Z d and that the final size of its outbreak is n; when n is fixed the function L n (p) can be regarded as the likelihood function for the percolation probability parameter p.) Letp n be the maximum likelihood estimate for p derived from L n (p). For any p define also L n (p) = L n (p)/L n (p c ).
Theorem III.1 The sequence of maximum likelihood estimatesp n for p converges to the critical probability p c .
The proof of this theorem is based on the following lemma, together with a generally accepted regularity condition.
Lemma III.2 For any p ∈ (0, 1) different from p c the following holds: Moreover, this convergence is uniform for any closed interval which does not contain p c .
The proof of the lemma and the theorem can be found in Appendix A 2.
As mentioned above, it is reasonable to expect the sequence of the maximum likelihood estimatesp n to be monotonic (see Figure 3). This observation gives rise to the following conjecture.
Conjecture III. 3 The sequence of {p n } converges to p c monotonically from the left.
Furthermore, we formulate the corresponding companion conjecture for the sequence of posterior distribution, using the notion of a delta sequence [20].
n=1 is a delta sequence which generates the delta function δ(p − p c ).
Thus, we believe that the limiting posterior distribution of the percolation parameter is a one-point mass distribution at p = p c , or the Dirac delta function δ(p − p c ).

Combinatorial characterisation of large percolation clusters on L d
The theoretical results obtained and conjectured previously for inference under scenario S2 can be used to derive their combinatorial analogues regarding the relative number of realisations of the process with the cluster size n. Under scenario S2 the posterior distribution π(p | |C| = n) can be seen as a mixture of beta distributions, as in (1). As in Section III B, let q(s, k, l) be the number of all possible graphs G which may correspond to an open cluster C (the vertex set of G and C coincide), and such that k is the number of edges in G, s is the number of edges in the saturation of G, and l is the number of edges of Π between the surface and the frontier of G. Conjecture III.4 would suggest that the number q(s, k, l) of graphs G corresponding to open clusters C which could emerge as a result of the percolation process with parameter p c and for which it holds that is far greater than the number of all other graphs. This is so, since the sequence of beta distributions Beta(α n , β n ) is a delta sequence generating the delta function at p c if and only if α n /(α n + β n ) → p c . Thus, in percolation processes on L d the number of finite graphs corresponding to open clusters of size n (where n is large) that satisfy the condition largely exceeds the number of all other connected components on n nodes. In other words, a typical graph corresponding to an open cluster of a large size in percolation process on L d is a connected subgraph of L d characterised by (3). In particular, when d = 2: that is the number of open edges in G is approximately equal to the total number of closed edges and edges between G and its frontier.

Large percolation clusters as rare events
When n is large, the appearance of finite open clusters of size n is highly unlikely: the distribution of the cluster size (hypothetically) decays as n −1−1/δ , δ > 0, when p = p c , and the decay is exponential (sub-exponential) when p < p c (p > p c ). Large finite percolation clusters can therefore be viewed as rare events. Since the state space of the McMC proposed for inference on the percolation parameter p under scenario S2 and described in Algorithm 2 involves the set of all open clusters on n nodes, this algorithm can be readily used in order to obtain realisations of these rare events.

IV. BAYESIAN OPTIMAL DESIGNS
A. Utility based optimal designs

General description
The experimental design problem can be conveniently approached within the Bayesian framework [13]. Suppose we study a stochastic process for which we formulate a model M , characterised by a model parameter θ. The model M specifies a sampling density f d (y | θ) for the experimental outcome y obtained from observation of the process under experimental set-up d given the value of the model parameter θ. Our knowledge about θ is described by a prior distribution π(θ). Whenever the choice of set-up d is under our control there arises the question of selecting the optimal d under which one should observe the stochastic process. Such prescribed conditions are referred to as a design, and the optimal design is found under optimality criteria specifically tailored to the context and the purpose of the experiment.
By employing a utility function u(d, y, θ) one can specify the purpose of the experiment and measure the value of its outcome y accordingly. The methodology of posing and solving utility-based optimal design problems within the Bayesian paradigm became somewhat standard [14,21]. The design has to be chosen before performing the experiment and one may choose to maximise the expectation of the utility function u(d, y, θ) with respect to θ and y [14]: where Here D is the set of possible designs. The set of possible outcomes y of the experiment is denoted by Y . The experiment is defined by a model f d (y | θ), i.e. by the sampling density of y conditional on θ for a given design d.
When the purpose of the experiment is to infer the model parameter, a sensible choice to measure the utility of the experiment under design d might be one that represents the information gained on θ from the experiment. One of the standard choices is the Kullback-Leibler (KL) divergence D KL between the prior distribution of θ, π(θ), and its posterior distribution π(θ | y, d) obtained following the experiment yielding data y: Introduced by Kullback and Leibler [22] this information measure was justified by Lindley [23] and Bernardo [24] and is related to the Shannon information measure.
The KL divergence is a random variable, being a function of the random outcome y. In order to find the the design d max that is optimal with respect to this choice of utility measure one should maximise the expected value of the KL divergence over the set of observables Y as a function of the design variate d: this being invariant under a reparametrisation of the model in terms of θ.
When the experimental motivation consists in increasing one's knowledge on the model parameter the expected KL divergence U KL (d) coincides with (6) taking u(d, y, θ) = log π(θ | y,d) π(θ) . If, however, the purpose of the experiment is to inform someone who holds a different prior to our own and we wish to advise them on which design to use, using our own 'superior' knowledge of the system under study, then the expected KL divergence should be calculated in the form prescribed by (8) when the integration over the space of observables Y is to be carried out using one's superior knowledge. We refer to these contrasting experimental scenarios as progressive design and instructive design [29], respectively, and discuss them next in more detail.

Progressive design
Under this scenario there is an experimenter E who holds prior beliefs regarding the parameter, θ, in the form of a prior distribution π(θ) (perhaps obtained from earlier observation of the process) and whose goal is to design an optimal experiment in order to increase this knowledge. Using KL divergence as the measure of information gain, this experimenter should maximise the expected KL divergence (8): so that U KL (d) coincides with (6) where u(d, y, θ) = log π(θ | y) π(θ) . Notice, however, that the expected information gain U KL (d) can be written as follows: where the second term is the Shannon entropy of π(θ), which does not depend on d. Thus, in order to obtain the solution d max in this case, it suffices to maximise the first term only.

Instructive design
In contrast to the progressive design scenario, in the instructive case there is an experimenter E , holding a prior π(θ), and a better (or differently) informed trainer (or instructor) T whose knowledge about the model parameter is summarised in a distribution π * (θ). (Of particular interest is a special case where the instructor knows the true value of θ.) The aim here is to maximise the change in experimenter's belief from π(θ) to π(θ | y) by designing an experiment using the superior knowledge π * (θ).
This latter optimisation problem can be formulated as follows: where, as before, π(θ | y) is derived as in (7), and f * (y) is the instructor's prior predictive distribution for y: In particular, if the instructor T knows the exact value of θ, θ * , and hence π * (θ) is the Dirac function δ(θ − θ * ), then f * (y) = f d (y | θ * ), so that Note that the 'instructive' expected utility U * KL is generally not symmetric with respect to the exchange of π(θ) and π * (θ): if the experimenter E and the instructor T change their roles, then the corresponding optimal instructive designs are different, unless their knowledge about θ is identical, that is π(θ) = π * (θ) almost everywhere.

B. Simulation-based evaluation of the expected utility
Generally, the solution to the optimal design problems (5-6) and (9-10) cannot be obtained analytically. This is mainly due to the following three problems. First, the design space D can be complicated, with many design variables, some of them having a continuous range of values. Second, even if the design space has a simple structure, the utility function may not be simple to evaluate, so the expected utility U (d) cannot be obtained explicitly. Finally, for incomplete observations of highly non-linear stochastic processes such as epidemic models, the likelihood is not usually available in a closed form, and this results in computationally intensive evaluations or sampling methods.
A review of analytical and approximate numerical solutions to Bayesian optimal design problems for traditional experimental design involving linear and non-linear models can be found in [13,15].
Müller [14] reviews simulation-based methods for optimal design problems where the expected utility U (d) is evaluated by Monte-Carlo simulation. In its simplest form an estimate U of U for any given design d in the progressive case is as follows: where {(θ i , y i ), i = 1, . . . , M } is a Monte-Carlo sample generated values: The expected utility U may, in particular, be based on the KL divergence. Analogously, the expected utility U * KL under the instructive scenario, when the instructor knows the true value of the model parameter, θ * , can be evaluated using the following scheme: where y * i ∼ f (y | θ * , d), i = 1, . . . , M , and KL(y i , d) is an estimate of the KL divergence D KL {π(θ | y * , d) π(θ)} that can be obtained via numerical integration of log π(θ | y,d) π(θ) with respect to the posterior π(θ | y, d).
The former function, in turn, may need to be evaluated through simulation methods-perhaps using an McMC scheme.
Evaluation of a continuous expected utility surface U (d) by (12) or (14) can be achieved by computing its values on a discretised grid of points and further smoothing of the obtained set of values in order to approximate the expected utility landscape. This, however, may be problematic when d ∈ D is a multidimensional design parameter. An alternative method is provided by the augmented probability simulation approach which is studied in [25,26], and reviewed in [14].
The augmented probability simulation approach assumes that u(d, θ, y) is a non-negative bounded function. (This condition, although not always automatically satisfied, can be easily achieved by correspondingly modifying the utility function.) An artificial density, proportional to u(d, θ, y)f d (y | θ)π(θ) can then be defined for (d, θ, y) jointly [14]: so that its marginal in d is proportional to the expected utility U (d). Sampling from h(·, ·, ·) can be used to obtain a sample from its marginal using the Metropolis-Hastings McMC scheme described in [14]. When using the augmented probability modelling approach one approximates the optimal design by an empirical mode of the marginal density of d under the artificial distribution h, that is by the mode of histogram formed from the first components of the set of samples obtained. Notice that the evaluation of u(d, θ, y) at each sampling step may in general itself involve an McMC sampling from the posterior π(θ | y, d), so that the approach is potentially very expensive computationally.

C. Inner-outer plots
We now apply the utility-based Bayesian approach introduced above in order to design optimal experiments for our SIR epidemic model with incomplete observations of the kind described by scenario S1. We do so by identifying a (finite) class of designs (node configurations) which we call inner-outer plots. These are plots of bounded size which are obtained by removing some nodes from the underlying grid. We refer to this process as sparsification of the grid.
A typical example of an inner-outer plot in Z 2 is given in Figure 4. The inner-outer plot depicted consists of two parts-inner and outer regions (hence the name of the design). Some sites are removed from the outer region of the plot in a regular manner, while all sites of the inner region are preserved. Our motivation in considering sparsified grids is to show that intuitively obvious designs, such as completely dense grids, can be improved. The justification of the grid sparsification is that it impedes the ability of the epidemic to spread, which makes possible better inference when p is large, while the more dense inner part of the design works better when p is N (m, r), r > 0, (16) where C (d) N is defined as follows [30]: We call any plot that can be obtained by translating the plot Π (d) 0 (m, r) in Z d an inner-outer (m, r)-plot, or simply an inner-outer plot, and denote it by Π (d) (m, r).
The total number of nodes contained in an (m, r)-plot can be calculated by subtracting the total number of the nodes removed from the outer plot (see (16)) as follows: where, as before, The inner-outer (m, r)-plot depicted in Figure 4 is from Z 2 . The size of the inner plot is m × m, where m = 9, and there are r = 3 'circles' (with respect to the metric · ∞ ) in the outer plot from which every second node is removed. The size of the bounding box is N × N , where N = m + 4r = 21. The total number of nodes T that this node configuration contains can be calculated using (17) and is equal to 357.

D. Inference for SIR epidemics on inner-outer plots
We consider an SIR epidemic with constant infectious times which arises at the central site of an inner-outer plot Π = Π (d) (m, r) from Z d and evolves on that plot according to the nearest-neighbour interaction rule, that is using edges from L d with endpoints from Π, subject to being restricted by the plot boundary. We consider this model in the context of scenario S1, that is when the only information available about the outcome of the epidemic is its site configuration C. As before, this model is equivalent to that of percolation on Π (d) (m, r).
We note that Algorithm 1 and Algorithm 2 can be readily used for making inference on the percolation probability parameter p given the sites of an open cluster C, or merely its size in the case of Algorithm 2, generated from a percolation process on any locally finite graph [31]. Figure 6 shows the plot of a simulated open percolation cluster C obtained on the inner-outer plot Π (2) (13, 2) with p = 0.9 (left plot) and an estimate of the likelihood formed from a histogram estimate of the posterior density of p obtained using the Markov chain described in Algorithm 1 with C as the input data, with a uniform prior density for p. Figure 5 shows similar plots for a realisation of the percolation process on Π (2) (23, 4) with p = 0.86.

E. Optimal design problem and design space
We consider now an optimal design problem for the percolation process on node sets of limited size. We adopt the utility-based Bayesian approach formulated in Section IV A. The choice of the design space D can be made in a number of ways, possibly reflecting such restrictions as limitations on the number of experimental units or the dimensions of the experimental plot. In the context of inner-outer plots the former restriction would mean  that T from (17) is bounded, whereas the latter condition is equivalent to bounding the quantity N (m, r). A combination of these conditions or some other information can be also taken into account when identifying the design space. One advantage of grid-based plot designs is that they can greatly reduce the dimensionality of the design space compared with designs that allow units to be placed at arbitrary locations.
We assume now that N is odd and fixed and that the design space has the form For example, if N = 19 then, as can be derived easily from (18), the design space consists of the following designs: The set of observables Y is the set of all connected components on d ∈ D containing the central node.

V. PRACTICAL IMPLEMENTATION
Given a finite design space D it is fairly straightforward to solve the optimal design problem for the percola-tion model on inner-outer plots using the tools developed above. We now make some comments about solving the problem under each of the two design scenarios.

A. Progressive design: expected utility evaluation through augmented modelling
Since D is finite we choose to identify the design that maximises the expected utility function in the progressive case using augmented modelling as described in Section IV B. Recall that this is based on an artificial distribution h(d, p, y) ∝ u(d, p, y)f d (y | p)π(p) from which samples are generated using a Metropolis-Hastings sampler. The optimal design d * is identified then as a value of d at which the marginal of h is maximised.

B. Instructive design: Monte Carlo evaluation of the expected utility
We treat the 'instructive' case differently from that of the 'progressive' one because of the form of the expected utility under this scenario. Recall that in the 'instructive' case whenever the instructor knows the true value p * of the model parameter p one can write the expected util-ity based on the Kullback-Leibler divergence as follows (compare with (11)): where f d (y | p * ) is the likelihood function evaluated at the true value of the model parameter p * given the open cluster y at the centre of the inner-outer plot d.
It is clear that to evaluate the expected utility U * KL (d) via standard Monte Carlo simulation, then a Markov chain has to be run to evaluate the posterior π(p | y, d) each time we sample for a new observation (an open cluster) y and also the potentially time-consuming integration has to be done with respect to the model parameter p. This integration can be implemented in the following way. Since we can sample from the posterior π(p | y, d) via McMC (Algorithm 1) for any given open cluster y, we do so and then fit the beta distribution (or some other distribution) to the McMC sample obtained in order to perform integration numerically in a more efficient, albeit, approximate manner.
Thus, the expected utility evaluation scheme for innerouter plots in the 'instructive' case and scenario S1 can be described as follows.
For each inner-outer plot d ∈ D we do the following: (ii) perform M McMC's in order to obtain M independent samples for the posterior distribution π(p | y, d); (iii) fit a Beta distribution to each of the samples obtained; refer to the fitted distributions as π(p | y i , d); (iv) evaluate numerically the integrals (v) estimate the expected utility: It is natural to ask how well the true posteriors can be fitted by a Beta distribution (step (iii) of the above scheme). Experience suggests that such a fit never affects the outcome of the analysis on the qualitative level unless the prior distribution has more than one local mode or its support is smaller than the entire interval (0, 1). Recall that the purpose of this step is to make evaluation of the integrals I i easier and faster, and hence the family of Beta distributions is only one possible choice. For instance, if the prior distribution π(p) is uniform on (0, 1), then the integrals I i represent the entropies of the fitted distributions. In the case when the fitting is done by beta distributions, each of these integrals can be quickly calculated using the following analytical formula for the entropy of the beta distribution Beta(α, β): where ψ is the digamma function, ψ(z) = Γ ′ (x)/Γ(z).
Other choices of prior distribution may suggest alternative and more appropriate fitting distributions to facilitate the calculation of the integrals I i , i = 1, . . . , M .

C. Example
In our example we consider all inner-outer plots in L 2 whose sizes do not exceed N = 11. There are only three such plots: Π (2) (3, 2), Π (2) (7, 1), and Π (2) (11,0). For ease of reference we mark them A, B, and C respectively (as depicted in Figure 7). Thus, the design space D = {A, B, C} consists of three designs, among which A is the most sparsified plot whereas no nodes are removed from C at all. Figure 8 represents graphically the results of the comparison of designs from D under both 'progressive' and 'instructive' scenarios when the prior distribution π(p) is uniform on the interval (0, 1). The left panel of the figure corresponds to the former scenario and depicts a histogram of a sample corresponding to the marginal of the artificial augmenting distribution h(d, p, y) ∝ u(d, p, y)f (y | p, d)π(p) in d ∈ D. The right panel corresponds to the latter scenario and shows the Monte Carlo estimated values of the expected utilities and 95% credibility intervals for each of the three considered designs (M = 1500, see (12) in Section IV B) assuming that the instructor knows p precisely so that π * (p) = δ(p − 0.9).
The plots from Figure 8 indicate that the solutions to the optimal design problem under the two scenarios are different from each other. The 'moderately sparsified' plot B maximises the expected utility in the progressive case, that is in the case when the experiment is designed by a single experimenter. If, however, it is the instructor who knows the true value of the model parameter (p = 0.9) and wishes to choose the best inner-outer plot from the set D for a more uninformed experimenter to use (instructive scenario), then the optimal plot is the 'mostly sparsified' inner-outer plot A. Notably, the 'mostly dense' plot C would be the worst choice in the instructive case, whereas it outperforms the 'mostly sparsified' plot A but is worse than the 'moderately sparsified' plot B in the progressive case. These results demonstrate the moderately counter-intuitive finding that the optimal design is not necessarily that with the maximal number of experimental units and that a 'less-is-more' principle may be in operation. On another level, this finding may not be entirely unpredictable in the instructive case. The optimal design's combination of a dense centre and sparsified outer region may help to increase the probability that the experimental epidemic succeeds in establishing itself but does not lead to saturation of the lattice before it dies out. This in turn may represent an outcome that yields much information on the unknown p.

VI. CONCLUDING COMMENTS
In this paper we have considered a number of question relating to the observation of percolation processes on finite lattices. We have demonstrated the value of McMC techniques for parameter estimation under different sampling regimes and shown how Bayesian experimental design techniques can be used to identify the optimal lattices on which to carry out experiments.
We introduced inner-outer plot designs by sparsifying the underlying grid. Allowing one to control excessive connectivity in outer regions of experimental plots, these designs help to separate different values of the model parameter p when p ≫ p c , while giving the epidemic a chance to establish in inner regions when p ≪ p c . Inner-outer plots can be further generalised by considering more than two levels of grid sparsification, according to the principle used in our work: the farther the region is from the centre of the plot the more intensive sparsification it receives. To an extent, this process can be viewed as switching to similarly shaped lattices with percolation parameter taking its values from the set {p, p 2 , p 3 , . . .} in decreasing fashion.
Although we have worked largely within framework of percolation theory, the connection between percolation and spatio-temporal epidemic models with nearestneighbour interactions mean that the results may have significant implications for and application to the study of epidemic systems.
Algorithm 1 Markov Chain Monte Carlo: scenario S1 Require: an open cluster C; 1: take an initial value p0 arbitrary from (0, 1); 2: t := 0 Xt := (pt, Satur C); 3: repeat 4: Gibbs sampler steps: 5: p2t+1 ∼ Beta (e(G2t) + 1, e(Satur C) − e(G2t) + wC + 1); 6: X2t+1 := (p2t+1, G2t); 7: Metropolis-Hastings sampler steps: 8: choose an edge e uniformly at random from Satur C; 9: if e ∈ G2t then 10: U ∼ Uniform[0, 1]; 11: if U ≤ min(1, this algorithm. Assume that the current graph within the Metropolis-Hastings step is G and a graphG is proposed, the latter being obtained from the former by deleting a vertex u with all edges adjoining it and inserting a vertex v, including each potential edge to v independently with probability p, this probability having been drawn in the preceding Gibbs step. We assume that at least one such edge is inserted and denote the number of deleted and added edges by d(u) and d(v) respectively. Then, and similarly, Clearly, so that the acceptance probability at the Metropolis-Algorithm 2 Markov Chain Monte Carlo: scenario S2 Require: the value of n. 1: take a value p0 arbitrary from (0, 1) and a graph G0 arbitrary from Gn; 2: t := 0 Xt := (pt, Gt); 3: repeat 4: move_on:=1; 5: Gibbs sampler steps: 6: p2t+1 ∼ Beta (e(G2t) + 1, e(Satur G2t) − e(G2t) + w(G2t) + 1); 7: X2t+1 := (p2t + 1, G2t); 8: Metropolis-Hastings sampler steps: 9: choose a vertex u uniformly at random from G2t+1 and choose a vertex v uniformly at random from ΓG 2t+1 . Derive a graphG from G2t+1 by deleting all edges which adjoin u (in G2t+1) and adding the edges that connect v with vertices of the graph G2t+1 \ {u} in Π, each independently with probability p2t+1 (conditioning on the event that at least one edge is added). 10: ifG is disconnected then 11: X2t+2 := (p2t+1, G2t+1); move_on:=0 12: if move on then 13: U ∼ Uniform(0, 1); 14:d(v)  Hastings step, α, is as follows: where, as it was introduced in the description of Algorithm 2, We claim that the constructed chain is irreducible. In terms of graph theory this means that the chain can proceed from each connected graph on n vertices, including the origin, to any other connected graph including the origin of the same size on the underlying lattice. We show the irreducibility of the proposed McMC by constructing a sequence of steps in which any graph of G n is transformed to a so called line-skeleton graph on n vertices. By such a graph we mean any tree (a graph with no cycles) containing the origin and having n vertices so that only two of these vertices have degree one. Select one such line-skeleton and denote it by S.
Consider a graph G from G n and denote the length of the shortest open path from x ∈ G to S by δ(x, S) (chemical distance). Each vertex x from G receives a well defined finite weight δ(x, S), since the graph G is connected and finite. By using the description of our Markov chain we can delete any vertex from our current graph for which δ(x, S) is maximal and add to this graph a vertex from the chosen line skeleton S without making the graph disconnected or containing cycles until the maximum value of d(x, S) is zero-in this case all vertices are forming the line skeleton. Since this procedure can be reversed it follows that G n in the described Markov chain is indeed a communicating class, and hence the chain is irreducible.

Lemma and Theorem
Proof. (Lemma III.2) We consider the following two cases [19]: (i) Subcritical case p < p c . In this case the cluster size distribution L n (p) decays exponentially, i.e.
However, it is also known that at p = p c the cluster size distribution, while almost surely finite, nevertheless has an infinite mean. This implies that, for any λ such that 0 < λ < λ p , we have lim sup n→∞ e λn L n (p c ) = ∞, which when combined with (A3), gives lim inf n→∞ L n (p) = 0, as required. Uniformity of convergence on any closed interval contained within (0, p c ) follows since λ above may be chosen so that λ p is uniformly bounded away from λ for p belonging to such an interval.
The infinite mean of the cluster size at p = p c also implies that for any η such that 0 < η < η p , we have lim sup n→∞ e ηn (d−1)/d L n (p c ) = ∞, which, when combined with (A4), again gives lim inf n→∞ L n (p) = 0, as required. Uniformity of convergence on any closed interval contained within (p c , 1) follows as before, and, when taken with the earlier result, gives the uniformity assertion of the lemma.

Proof.
(Theorem III.1) We assume only that lim n→∞ L n (p) exists. Then, by Lemma III.2 this limit is equal to zero for all p = p c . (The existence of the limit would follow, for example, from the widely believed, but not rigorously proved, result that for p = p c the distribution of the size of the open cluster size containing the origin has a power-law tail-see Chapter 9 of [19].) Then the result of Lemma III.2 can be reformulated in ε-terms as follows: ∀ε > 0 ∃ N (ε, γ) > 0, such that L n (p) < εL n (p c ) ∀n > N (ε, γ) ∀p ∈ [α, p c −γ]∪[p c +γ, β), (A5) for any α and β, such that ∆ = [α, p c − γ] ∪ [p c + γ, β], ∆ ⊂ (0, 1).
Appendix B: Mixing properties of chains from the proposed algorithms Scenario S1 and Algorithm 1. Figure 9 displays a trace plot of updates in p for inference made for the site configurations from Figure 1 (Algorithm 1 was applied) and presented in Figure 2. The trace plot indicates that the mixing properties of the chain are rather satisfactory. The Matlab-based simulation used required 15 seconds on Intel(R) Core(TM)2 Duo CPU 2.26GHz to obtain a series of chain updates of the length 10 4 .
Experiments with larger configurations suggest that the burn-in time of the constructed chain may vary considerably with the size of the site configuration and its potential connectivity properties as well as on the choice of the initial graph G 0 (in our experiments we worked with fully saturated graphs and with connected components derived from those by random sparsification of edges). However, in practice, the chain described by Algorithm 1 typically converges rapidly: both the connected component updates and the percolation parameter updates can be efficiently realised, allowing one to  Figure 2 for the cluster C in Figure 1.
generate a suitable sample from the posterior distribution π(p | C) within practically useful timescales. Updates of the connected component, that is deletion and insertion of edges while preserving the connectedness of the underlying graph, can be further improved using dynamic graph algorithms [27]. Scenario S2 and Algorithm 2. Figure 10 shows histograms of samples from the distribution π(p | |C| = n) obtained using the proposed McMC for the scenario S2 for the cluster size values p = 10, 35, 50, 70 and corresponding trace plots.