# Prediction and Surveillance Sampling Assessment in Plant Nurseries and Fields
## 2 Method
### Implementation of Bayesian Kriging with INLA
Matern fields are the only stationary solutions to a stochastic partial differential equation (SPDE). The goal of the SPDE approach is to find a Gaussian Markov random field (GMRF) with a neighborhood structure and a precision sparse matrix \(Q\) that best represents the Matern field. Given this representation, inferences can be made using the GMRF found and its good properties of calculation. Basically, the SPDE approach uses a finite representation to define the Matern field as a linear combination of base functions, which are defined by a triangulation into domain \(D\). This triangulation divides \(D\) into a set of triangles that are not intercepted and united by at least one common edge or corner. First, the initial vertices of the triangles are placed at locations \(s_{1},\ldots,s_{n}\), and later, additional vertices are added in order to obtain a useful triangulation desired for spatial prediction. Considering the triangulation, the representation of the base function of a Matern field \(X(s)\) is given by  
\[X(s)=\sum_{i=1}^{n}\psi_{l}(s)w_{l} \tag{10}\]  
where \(n\) is the total number of vertices, \(\{\psi_{l}(s)\}\) are the basis functions and \(w_{l}\) are weights with a Gaussian distribution. The functions \(\{\psi_{l}(s)\}\) are selected as linear segments in each triangle, i.e., \(\psi_{l}(s)\) is 1 on vertex \(l\) and 0 on the other vertices. The height of each triangle (the space field value at each vertex of the triangle) is given by the weight \(w_{l}\), and the values inside the triangle are determined by linear interpolation. The key point of the SPDE approach is the finite representation (10) that establishes the link between Gaussian fields (GF) \(X(s)\) and the GMRF defined by Gaussian weights \(w_{l}\). These weights may be assigned a Markovian structure as shown in [9]. In particular, the precision matrix \(Q\) of GMRF is defined by the equation \(w_{l}\)\(\sim\)\(N(0,Q_{S}^{-1})\) as a function of \(\kappa^{2}\) for \(\alpha=1,2,\ldots\), \(\
u=0,1,2,\ldots\) and \(\alpha=\
u+1\). Under this perspective, for each vertex \(i=1,\cdots,n\), the full hierarchical model structure can be stated as follows:  
\[ \begin{array}{c}Y_{i}\sim\text{Ber}(\mathbb{g}_{i})\\ \text{logit}(\pi_{i})=\beta_{0}+W_{i}\\ \pi(\beta_{0})\propto 1\\ W\sim N(0,\sigma_{W}^{2}H(\phi))\end{array}  \tag{11}\]  
In contrast to WinBUGS software [18], regarding assignation of the priors, the correlation function is not modeled directly. In this case, the numerical solution to the Gaussian field is made using an approximate weak solution to a stochastic partial differential equation (SPDE) as a Gaussian Markov random field [4; 15]. This solution requires defining two new parameters, \(\kappa\) and \(\tau\), which determine the range of the spatial effect and the total variance. More precisely, the range is approximated by the expression \(\phi\approx\sqrt{8}/\kappa\), while the variance is defined as \(\sigma_{W}^{2}=\frac{1}{4\pi\kappa^{2}\tau^{2}}\). \(\kappa\) and \(\tau\) have the following prior distributions:  
\[\begin{array}{c}2\text{log}\kappa\sim N(m_{\kappa},\sigma_{\kappa}^{2})\\ \text{log}\tau\  
reasonably according to the size of the region, while the mean \(m_{\tau}\) is chosen so that the variation of the field is 1.