Image Anal Stereol 2014;33:75-81 doi:10.5566/ias.v33.p75-81 Short Research Communication MODELS OF COVARIANCE FUNCTIONS OF GAUSSIAN RANDOM FIELDS ESCAPING FROM ISOTROPY, STATIONARITY AND NON NEGATIVITY PABLO GREGORIC,1, EMILIO PORCU2 AND JORGE MATEU3 1IMAC, Instituto de Matem´aticas y Aplicaciones de Castell´on, Departamento de Matem´aticas, Universitat Jaume I de Castell´on de la Plana, Spain; 2Universidad Federico on, ESTCE, Campus Riu Sec, E-12071 Castell´Santa Mar´ia, Departamento de Matem´aticas y aticas at Valparaiso, Chile; 3IMAC, Instituto de Matem´Aplicaciones de Castell´on, Departamento de Matem´aticas, Universitat Jaume I de Castell´on, ESTCE, Campus Riu Sec, E-12071 Castell´on de la Plana, Spain e-mail: gregori@uji.es, porcu@usm.cl, mateu@uji.es (Received October 25, 2013; revised December 20, 2013; accepted January 18, 2014) ABSTRACT This paper represents a survey of recent advances in modeling of space or space-time Gaussian Random Fields (GRF), tools of Geostatistics at hand for the understanding of special cases of noise in image analysis. They can be used when stationarity or isotropy are unrealistic assumptions, or even when negative covariance between some couples of locations are evident. We show some strategies in order to escape from these restrictions, on the basis of rich classes of well known stationary or isotropic non negative covariance models, and through suitable operations, like linear combinations, generalized means, or with particular Fourier transforms. Keywords: anisotropy, covariance model, Gaussian random .eld, non negativity, non stationarity. INTRODUCTION A typical assumption in image analysis is that of the white noise: a Gaussian error distorts the true value of each pixel. Other reasonable approaches admit that perturbations affecting nearby pixels can be correlated in distinct manners (Winkler, 2006), and our work deals with processes within this framework. Although digitized images form lattices of pixels, one can analyze the behavior of randomness in large images, with an acceptable degree of approximation, by considering images as continuous domains, hence having at hand the tool of random .eld models, from (space or space-time) Geostatistics (Cressie and Wikle, 2011). Once the true value of the image is subtracted, the noise can be reasonably modeled as a zero mean Gaussian random .eld (GRF), which can be completely determined by the so called covariance function. A good estimation of the covariance function leads to the understanding of the noising process affecting the images taken by a particular device or technology. This step of estimation from sampled data is a rather complicate task, since one needs to have a rather dense (in the colloquial sense) set of observations in the whole domain. This is unfeasible in most of the cases, and the reason why additional simplifying assumptions are imposed on the GRF. As a consequence, such classes of GRFs may not represent well the behavior of processes in some .elds of application. One of these assumptions is weak stationarity (Cressie and Wikle, 2011), which stands for GRF such that the covariance among two locations depends only on the (vector) difference of locations. In this case, the covariance function can be written as a function of one (lag) vector, and all couples of locations with difference close to the lag serve to estimate the covariance function at that jump vector. Another assumption is isotropy (or radial symmetry) (Schoemberg, 1938; Daley and Porcu, 2013): covariance among two locations depends solely on their Euclidean distance. In this case, the covariance function can be written as a function of one non negative number, and all couples of locations showing similar distance are used to estimate the covariance function at that distance. There are several natural phenomena where the spatial direction plays a dominant role, and isotropy does not allow to take this fact into account. Throughout this paper, {Z(s) : s . Rd } shall denote the zero mean GRF, and C its covariance function, de.ned as C(s1,s2) := cov(Z(s1),Z(s2)). Any function K : Rd × Rd › R is the covariance function of some GRF if and only if it ful.lls the positive de.niteness condition, i.e., if for every set of locations {si}ni=1 . Rd of any size n . N, and every set GREGORI P ET AL: Covariances of GRF escaping from isotropy, stationarity and non negativity of coef.cients {ai}ni=1 . Rthe inequality nn .. aiK(si,sj)aj . 0 (1) i=1 j=1 holds, i.e., the matrix (K(si, sj))ni, j=1 is positive de.nite for every choice of {si}n i=1. Since GRFs are completely determined by their .rst and second order moments (Stein, 1999), parametric models of zero mean GRFs can be given by parametric families of positive de.nite functions. Under the assumption of stationarity for the GRF, we have that C(s1,s2)= C(s1 - s2,0) for all s1,s2. In this case it is usual to make the abuse of notation C(h) := C(s1,s2) whenever s1 - s2 = h. If isotropy is also imposed, then C(s1,s2)= C(s3,s4) if Is1 -s2I = Is3 -s4I. The abuse of notation for this case is C(h) := C(s1,s2) whenever Is1 - s2I = h. Literature provides several .exible and rich parametric models of isotropic covariance functions for GRF, such as the Gaussian (or squared exponential) model C(h; .)= .2 exp(-.h2) , (2) (. > 0), the Mat´ ern, 1960) ern class (Stein, 1999; Mat´ C(h;.,.) := .2(.h). K. (.h) , (3) (. > 0,. > 0 and K. a modi.ed Bessel function of second kind of order .), which includes the exponential model (. = 1/2, for unidimensional Ornstein–Uhlenbeck processes), the autoregressive of the .rst order model (. = 3/2) and even the Gaussian model (when . › .), several models of piecewise polynomial functions with compact support (Wendland, 2005), etc. All these families contain exclusively non negative covariance functions. It is obvious that previous models should not be used “as­is” in applications where negative covariances are natural (environmental processes, turbulences, etc.). In a particular real life process under analysis, any of these properties (or assumptions or features) might be discredited either by the nature of the process or due to suspicion after some graphical analysis of the dataset. In order to escape from isotropy, stationarity and non negativity, the common philosophy is to use the mentioned parametric families as building blocks for the construction of more sophisticated models. In this paper we make a review on the attempts for obtaining classes of covariance (or correlation) functions, escaping from the typical restrictions using rather simple strategies, such as projections of the Euclidean space, generalized means, linear combinations with negative coef.cients, etc. of simple ingredients. Our review does not pretend to be exhaustive of the whole literature, but instead focuses on the effort made in a line of research by this group of authors (Porcu et al., 2007; Gregori et al., 2008; Porcu et al., 2008; 2009). These new classes are hence suitable for .lling in some gaps in the modeling task. In the following sections we describe the general setting, a brief summary of other authors’ approaches, and our contribution in the direction of relaxing each of the considered assumptions. We end up by justifying the past and future efforts that shall be devoted into this line of research. ESCAPING FROM ISOTROPY Isotropic covariance functions can be written in terms of functions . de.ned on [0,.), such as C(h) := .(IhI) or, equivalently, C(h) := .(IhI2). They are sometimes referred to as radial functions. Matheron (1965) proposes to start with an isotropic covariance function C0, i.e., C0(h) := .(IhI2), and use a positive semide.nite quadratic form Q, given by Q(h) := hTMh (being M the matrix de.ning Q), in order to de.ne C := . . Q and escape from isotropy. The dimension d, up to which one can lift such a function . to a valid isotropic covariance function C in Rd , depends highly on the behavior of .. When using Euclidean norms, isotropic covariance functions, for each dimension, are perfectly determined by the classical result of Schoemberg (1938), written in terms of positive de.nite functions, but rephrased here in terms of covariance functions. Theorem 1 (Schoemberg, 1938). C(h) := .(IhI2) is a (isotropic by construction) covariance function in Rd ... – ... for a speci.c d if and only if there exists a non negative bounded measure F such n . that .(t)= 0 .d(ts)dF(s), where .d (t) := .(d/2)(2/t)(d-2)/2J(d-2)/2(t) (being J(d-2)/2 the Bessel function of the .rst kind). In other words, . it is the Fourier–Bessel or Hankel transform of such an F. – ... for all d . N if and only if there exists a non negative bounded measure F such that .(t)= n. exp(-ts)dF(s) (or equivalently if and only if 0 . is a completely monotone function). In other words, . it is the Laplace transform of such an F. Image Anal Stereol 2014;33:75-81 Functions . de.ning valid isotropic covariance functions in Rd are said to form a class denoted by .d (Gneiting, 1999; Daley and Porcu, 2013). It is easy to see that these classes are nested, .1 . .2 . ··· . .., where .. stands for the class of completely monotone functions. Matheron (1965) and Gneiting (2002) de.ne and use operators allowing to obtain new isotropic covariance functions in one of theses classes by transforming existing isotropic covariance functions lying in another of these classes. Indeed: Proposition 2 (Gneiting, 2002). The operators mont´ ee and descente, de.ned as [. › I.] and [. › D.] with: n . u.(u) du .'(t) t I.(t) := n . , D.(t) := (4) 0 u.(u) du t.'(0) .t > 0 and D.(0) := 1, satisfy the following properties: – If . . .d , d . 3 and the mapping [u › u.(u)] is integrable then I. . .d-2. – If . . .d and ..''(0) then D. . .d+2. With these tools, Porcu et al. (2007) proposed to project the whole Euclidean domain into arbitrary Cartesian projections in order to abandon isotropy, but preserving isotropic blocks. For that aim, let us consider Rd = Rd1 × ··· × Rdk with vector d . (d1,...,dk) determining the projected dimensions – hence d1 + ··· + dk = d – and denoting Rd : h . (h1,...,hk) with hi . Rdi . Now classes .d are generalized to: .d := {. : [0, .)k › R : C(h)= .(Ih1I,...,IhkI) is a cov. fct. in Rd }, (5) and their characterization follows easily: . . .d iff . . k .(s1,...,sk)=···..di (siui)dF(u1,...,uk) 00 i=1 (6) for some non negative bounded measure F. Now one has the possibility of applying mont´ee and descente operators by partial differentiation or integration at a particular projection, getting a wider range of possibilities. If . =(.1,...,.k) .{0,1}k, the . -mont´ees ee and . -descente are compositions of mont´and descentes in the proper dimensions (at indices j where .j = 1), and the properties of mapping between classes .d are generalized. Concretely, I. . =(I.1 . 1 I.2 . ···. I.k ). and D. . =(D.1 . D.2 .···. D.kk )., 2 k 12 where I1 and D1 ee and are, respectively, the mont´ jjdescente operator with respect to the j-th variable, and I0 j and D0 j are both the identity map, for j = 1,2,...,k. Proposition 3 (Porcu et al., 2007). Operators . ­mont´ ee and . -descente satisfy: – If . . .d and . .{0,1}k with 2 + 2. . d and .I.1 . I.2 .···. I.k .(0) then I. . . .d-2. . 12 k – If . . .d and ... ' .(0) .. '. 2. (in the sense of the limit) then D. . . .d+2. . The application of this technical result is fruitful, as it allows for the de.nition of space-time covariance functions isotropic only in the space block, but not in the whole dimension domain. Example 4. If we take f any non negative integrable function in R2 we get a non negative bounded measure n F by de.nite integration, i.e., F(·)= (·) f dudv. Now we can use the kernel .1(s)= cos(s) in order to de.ne . . .(s,t) :=f (u,v)cos(su)cos(tv)dudv , (7) 00 and get . . .(1,1) (see Eq. 6), i.e., C(h1,h2) := .(|h1|,|h2|) is a valid non isotropic covariance function in R2. If not of direct interest, it helps to de.ne a space-time covariance function which is isotropic in (tridimensional) space and symmetric in time, but non isotropic, seen as a function of four dimensions. Indeed: C((h1,h2,h3),t) := D(1,0).(I(h1,h2,h3)I,|t|) , (8) where . is the function of Eq. 7, is such a covariance function. ESCAPING FROM STATIONARITY A standard way of building non stationary GRF is through a deformation of the domain. If C is a stationary covariance function in Rn, and u : Rd › Rn is any map, then C . u is in general non stationary in Rd . See an application of this approach in Sampson and Guttorp (1992). Paciorek and Schervish (2004) build a non stationary covariance in Rd using an isotropic stationary covariance (given by a completely monotone function), and a d × d matrix which allows for spatially adaptive parameters. Higdon et al. (1999) build non stationary RF using spatially varying kernels Ks : Rd › R (where s is the spatial location parameter) and white noise RFs ., as n Z(s) := RdKs(u).(u)du. In this case the covariance n function is C(s1,s2) := RdKs1 (u)Ks2 (u)du, and the GREGORI P ET AL: Covariances of GRF escaping from isotropy, stationarity and non negativity process can be chosen to be non Gaussian (Guttorp and Schmidt, 2013). In a similar way, Pintore and Holmes (2004) have de.ned a procedure based on the spatial adaption of spectral densities of parametric families of stationary GRFs in order to get new non stationary GRFs. The following classical statement is the key for another strategy in the construction of classes of covariance functions. Theorem 5 (Bochner, 1933). A continuous function C in Rd is a (stationary) covariance function if and only if there exists a .nite non negative measure F such that C(h)= eih·. dF(. ) . (9) Rd If C is integrable, the measure F is absolutely continuous with respect to the Lebesgue measure, and there exists a non negative function f with -i. f (. )= 1e·hC(h)dh , (10) (2.)d Rd which is called the spectral density of C and, whenever it exists, one has C(h)= eih·. f (. )d. . (11) Rd New (parametric) classes of GRFs can be given through (parametric) classes of non negative functions. Fourier transforms of these classes are valid covariance functions, and the limit to this methodology is the availability of their closed form. For well known families of covariance functions, such as Gaussian or Mat´ern class, the spectral density has a nice closed form. They become perfect candidates for the de.nition of new classes under this methodology, enlarging the family of models of GRFs at hand. As it was said before, new families of covariance functions were proposed by Pintore and Holmes (2004), based on a parametric family of spectral densities and using the brilliant idea of adapting the parameter of the spectral density to each different location. A convenient inverse (similar to Fourier) transform leads the way to non stationarity. Proposition 6 (Pintore and Holmes, 2004). Let { f. }... be a parametric family of spectral densities (of a parametric model of stationary GRF). Let [s › .(s)] be any mapping of locations to parameters. Then the construction i. ·(s1-s2) f 1/2 (. ) f 1/2 C(s1,s2) := e(. )d. .(s1).(s2) Rd (12) de.nes a valid (non stationary) covariance function if and only if f.(s) is absolutely integrable for all s. An application of this result to the classical Mat´ern family with parameter (.,.), and adapting [s › .(s)] yields a closed form for a non stationary GRF, with .(s) a spatially varying smoothness parameter. Indeed, they obtain a model of the type K(s1,s2)= C(Is1 - s2I,.,.(s1,s2)) for C as de.ned through Eq. 3. This process is dubbed spatially adapted Mat´ern family of parameter (.,.(·)), and it can be preferred to .t some speci.c data whenever the Mat´ ern model, with .xed smoothing parameter, seems inadequate. In Porcu et al. (2009), the authors have analyzed f 1/2 the role of the geometric mean (. ) f 1/2 (. ) .(s1).(s2) appearing in Eq. 12, and have found that a general .­mean (called Archimedean mean therein), i.e., A. ( f1, f2)(. ) := .(0.5.-1( f1(. )) + 0.5.-1( f2(. ))) (13) for . completely monotone, extends the previous result and provides new examples, such as the closed form of another spatially adapted Mat´ern family, this time with parameter (.(·), .), adapting the scale parameter of the stationary Mat´ ern class to have a different behavior along the domain. This procedure enlarges, in a similar way, the family of available models for non stationary processes, giving more chances to understand other types of covariance structures. Namely, the authors propose a model of the type K(s1,s2)= C(Is1 - s2I,.(s1,s2),.) for C as de.ned through Eq. 3. ESCAPING FROM NON NEGATIVITY As it has been mentioned in the introduction, most of the classical examples of covariance functions are parametric families of non negative functions. However, a classical example of oscillatory covariance function is C(h) := c(.IhI)-. J. (.IhI), based on the Bessel function of the .rst kind J. , which is valid for . . (d - 2)/2,. > 0 (Yaglom, 1987). Other families of covariance functions attaining negative values have been obtained in Gneiting (2002), starting from families of non negative compactly supported covariance functions of Wendland (1995), and applying the turning bands operator of Matheron (1973) (see Fig. 1). 78 Image Anal Stereol 2014;33:75-81 Fig. 1. Gneiting family obtained by using the turning bands operator on a Wendland’s class. Usual ways of creation of new covariance functions from existing ones are, for instance, linear combinations with non negative coef.cients – or mixtures – pointwise products, and pointwise limits. But in all these cases, if the ingredients are non negative so the result is. The intuitive way that we have explored is the consideration of linear combinations of the form C. = (1 - . )C0 + .C1 with C0 and C1 non negative. For . . [0,1] we have a convex combination, and we lie on the cone of non negative covariance functions. Outside [0,1], we intend to leave that cone, but stay inside the covariance functions class. We have been inspired by Ma (2005), where the author deals with linear combinations of speci.c families using a tool: the non negativity of the spectral density. We posed the general problem and got the following result. Proposition 7 (Gregori et al., 2008). If C0 and C1 are continuous integrable covariance functions, then C. :=(1 - . )C0 + .C1 is a valid covariance function if and only if (1 - max(1,M))-1 . . . (1 - min(1,m))-1 (14) where m := inf. f1(.)/ f0(.), and M := sup. f1(.)/ f0(.), and f0,f1 are the spectral densities of C0,C1, respectively, . This result can be applied to families having a known analytic expression of the spectral density, and whose dependence on the frequency . comes through its norm (for instance, Gaussian and Mat´ern), and the optimization computation of m and M results feasible. Example 8. Fig. 2 shows particular linear combinations of two (non negative) Gaussian covariance functions which .nally attain negative values. Fig. 3 shows realizations of GRF with covariance which is linear combination of Gaussian covariances. Fig. 2. C. (u) :=(1 - .)e-0.5u2 + . e-u2 for . = 2 (cont), 1.75 (dashed) and 1.5 (dotted). A very similar result is found for the space-time model C. (h,t) :=(1 - .)C0s(h)C0t (u)+ .C1s(h)C1t (u). CONCLUSIONS As far as space or space-time sampling in large size lattices or continuous domains, such as in image analysis, becomes more and more feasible, more complex models, not restricted to be stationary or isotropic, or even to have only non negative values, can be used in applications. Besides, speci.c applications or datasets might suggest that stationarity, or isotropy should be rejected, and may give hints on which one of the mentioned models of this survey might be used in order to explain the data. We think that models free of such restrictions have chances of success in data explanation, and defend that their conception should be based mainly on the nature of the process, but also on simplicity, in the sense of: – simple operations – such as linear combinations, projections, integral transforms – and – acting on simple restricted models – such as isotropic, stationary, well behaved families – such as the Mat´ ern one. GREGORI P ET AL: Covariances of GRF escaping from isotropy, stationarity and non negativity Fig. 3. Realizations of zero-mean GRFs with covariance function (upper left) C0, (upper right) C1, (lower left) . 0.5C0 + 0.5C1 and (lower right) 2C0 -C1, where C0 and C1 are Gaussian covariance functions with . = 0.1 . and 0.2 respectively (see Eq. 2). Only the last case has a covariance with negative values, which is translated into a realization with negative correlation among nearby locations. ACKNOWLEDGEMENTS The authors acknowledge two anonymous referees, for helping at the improvement of the paper, and specially to Moreno Bevilacqua, for his diligent assistance. He is the coauthor of R package CompRandFld (Padoan and Bevilacqua, 2013), and has added the family of linear combinations of Gaussian covariance functions for making possible the simulations shown in Fig. 3 The second author acknowledges the support of Chilean Proyecto Fondecyt Regular number 1130647. First and third authors acknowledge the support of Spanish projects of Ministerio de Ciencia e Innovaci´ on MTM2010-14961 and Universitat Jaume I de Castell´ on P1-1B2012-52. The authors are also grateful to the organizing and scienti.c committees of the 11th European Congress of Stereology and Image Analysis 2013 (11th ECS) held at Kaiserslautern, where the topic of this paper was presented. REFERENCES Bochner S (1933). Monotone Funktionen, Stieltjes Integrale und Harmonische Analyse. Math Ann 108:378–410. Cressie N, Wikle C (2011). Statistics for spatio-temporal data. Hoboken: Wiley. Daley DJ, Porcu E (2013). Dimension walks through Schoenberg spectral measures. P Am Math Soc. In press. Gneiting T (1999). On the derivatives of radial positive de.nite functions. J Math Anal Appl 236:86–93. Gneiting T (2002). Compactly supported correlation functions. J Multivariate Anal 83:493–508. Gregori P, Porcu E, Mateu J, Sasvari Z (2008). On potentially negative space time covariances obtained as sum of products of marginal ones. Ann I Stat Math 60:865–82. Guttorp P, Schmidt AM (2013). Covariance structure of spatial and spatio-temporal processes. WIREs. Computational Stat 5:279–87. Image Anal Stereol 2014;33:75-81 Higdon DM, Swall J, Kern J (1999). Non-stationary spatial modeling. In: Bernardo JM, Berger JO, Dawid AP, Smith AFM, Eds. Bayesian Statistics 6. Proc 6th Valencia Int Meeting 761–8. New York: Oxford University Press. Ma C (2005). Linear combinations of space-time covariance functions and variograms. IEEE Trans Signal Proc 53:857–64. Matern´B (1960). Spatial variation. Meddelanden fr°an Statens Skogsforskningsinstitut, Alm¨anna F¨orlaget, Stockholm. 49(5):1-144. (2nd Ed (1986), Lect Not Stat 36. Berlin: Springer.) Matheron G (1965). Les variables r´ees egionalis´et leur estimation. Paris: Masson. Matheron G (1973). The intrinsic random functions and their applications. Adv Appl Prob 5:439–68. Paciorek C, Schervish MJ (2004). Nonstationary covariance functions for Gaussian process regression. In: Thrun S, Saul L, and Sch¨olkopf B, Eds. Advances in Neural Information Processing Systems 16. MIT Press. Padoan S, Bevilacqua M (2013). CompRandFld: Composite-likelihood based analysis of random .elds. R package, version 1.0.3. http://CRAN.R-project.org/ package=CompRandFld Pintore A, Holmes C (2004). Non-stationary covariance functions via spatially adaptive spectra. Tech Rep. http://www.stats.ox.ac.uk/~cholmes/Reports/ spectral tempering.pdf Porcu E, Gregori P, Mateu J (2007). La descente et la montee´´ etendues: the spatially d-anisotropic and the spatiotemporal case. Stoch Env Res Risk A 21(6):683– 93. Porcu E, Gregori P, Mateu J (2008). Recent advances to model anisotropic space-time data. Stat Method Appl 17(2):209–23. Porcu E, Gregori P, Mateu J (2009). Archimedean spectral densities for nonstationary space-time geostatistics. Stat Sinica 19(1):273–86. Sampson PD, Guttorp P (1992). Nonparametric estimation of nonstationary covariance structure. J Am Stat Assoc 87:108–19. Schoenberg IJ (1938). Metric spaces and positive de.nite functions. T Am Math Soc 44(3):522–36. Stein ML (1999). Interpolation of spatial data. New York: Springer. Wendland H (1995). Piecewise polynomial, positive de.nite and compactly supported radial functions of minimal degree. Adv Comput Math 4:389–96. Wendland H (2005). Scattered data approximation. Cambridge monographs on applied and computational mathematics. Cambridge University Press. Winkler G (2006). Image analysis, random .elds and Markov chain Monte Carlo methods: a mathematical introduction. Stochastic modelling and applied probability. New York: Springer. Yaglom AM (1987). Correlation theory of stationary and related random functions. Volume I: Basic results. Springer.