Conditional Fractional Gaussian Fields with the Package FieldSim

Abstract:

We propose an effective and fast method to simulate multidimensional conditional fractional Gaussian fields with the package FieldSim. Our method is valid not only for conditional simulations associated to fractional Brownian fields, but to any Gaussian field and on any (non regular) grid of points.

Cite PDF Tweet

Published

April 3, 2016

Received

Mar 6, 2015

Citation

Brouste, et al., 2016

Volume

Pages

8/1

38 - 47


1 Introduction

Rough phenomena arise in texture simulations for image processing or medical imaging, natural scenes simulations (clouds, mountains) and geophysical morphology modeling, financial mathematics, ethernet traffic, etc. Some are time-indexed, some others, like texture or natural scene simulations, should be indexed by subsets of the Euclidean spaces R2 or R3. Recent data (as the Cosmic Microwave Background or solar data) are even indexed by a manifold.

The fractional Brownian motion (fBm), introduced by (and developed by ) is nowadays widely used to model this roughness. Fractional Brownian motions have been extended in many directions: higher dimensions with fields, anisotropy, multifractionality, etc. This paper is devoted to a simulation method for conditional Gaussian fields. This could improve, in the future, natural scene simulations by fixing for instance the valleys.

The simulation of fractional Gaussian processes is not difficult in dimension one (see a review of ). Let us recall the numerical complexity of some classical methods: the Cholesky method has a complexity of O(N3) where N is the size of the simulated sample path. For specific stationary processes (on a regular grid) the Levinson’s algorithm has a complexity of O(N2logN) and the Wood and Chan algorithm (see ) a complexity of O(NlogN).

In higher dimensions, the Wood and Chan method has been extended to stationary increments fields with the Stein’s method ; the fractional Brownian field can therefore be simulated on a regular grid of the plane. For general Gaussian fields on a general discrete grid, the Cholesky method is costly and exact simulations are no longer tractable. Approximate methods have been intensively developed (midpoint, ; turning bands, ; truncated wavelet decomposition) but for specific fields. On manifolds, simulation procedures based on truncated series of eigenfunctions of the Laplace-Beltrami operator are discussed in .

Our approach, presented in , is based on a 2-steps method with an exact simulation step plus a refined fast step, that is an improvement of the midpoint method. It has been implemented in the FieldSim package . The fieldsim simulation method can be applied to general Gaussian processes on general simulation grids (regular and non regular) on Euclidean spaces and even on some manifolds (see Figure 1). It is worth mentioning that another package, RandomFields , allows the simulation of a large class of random fields such as Gaussian random fields, Poisson fields, binary fields, chi-square fields, t fields and max-stable fields (see ). In RandomFields, conditional random fields (which are the purpose of the present paper) are given for a wide range of spatial and spatio-temporal Gaussian random fields. Some of the default models of the FieldSim package cannot be simulated with the help of default models of the RandomFields package. Nevertheless, it is still possible to simulate them with the RMuser() and RFsimulate() commands of the RandomFields package. It may be noted that the FieldSim package does not allow for the simulation of more than the RandomFields package. FieldSim package is an alternative in which the underlying methods of simulation are generic.

graphic without alt textgraphic without alt text

Figure 1: On the left: fractional Brownian field (top-left), multifractional Brownian field (bottom-left), fractional Brownian sheet (top-right) and hyperbolic fractional Brownian field (bottom-right); on the right: fractional Brownian field on the sphere.

We propose here to adapt the FieldSim package to conditional simulations. Definitions and notation will be introduced in the following section with the “process” class, the setProcess procedure and the fieldsim procedure. The fieldsim procedure adapted to conditional Gaussian fields is described in the next section. Simulations with the package FieldSim are presented in the last section.

2 Notation and preliminaries

Fractional Gaussian fields

Let d be a positive integer and X()={X(M),MRd} be a real valued non stationary field with zero mean and second order moments. It is worth emphasizing that we consider in this paper the metric space Rd with the Euclidean norm but the method can be generalized to a smooth and complete Riemannian manifold equipped with its geodesic distance .

The covariance function R(,) is defined by: R(M1,M2)=cov(X(M1),X(M2)),M1,M2Rd. This function is nonnegative definite (n.n.d.). Conversely, for any n.n.d. function R(,), there exists an unique centered Gaussian field of second order structure given by R(,).

Different classical fractional Gaussian fields have been simulated to illustrate the FieldSim package in . In the sequel, M and M are two points of Rd and is the usual norm on Rd, d=1, 2. We can cite:

  1. The standard fractional Brownian fields are defined through their covariance function (e.g., ): R(M,M)=12(M2H+M2HMM2H), where the Hurst parameter H is real in (0,1).

  2. The standard multifractional Brownian fields are defined through their covariance function (see ):

    R(M,M)=α(M,M)(MH~(M,M)+MH~(M,M)MMH~(M,M)), where H~(M,M)=H(M)+H(M),α(M,M)=C(H(M)+H(M)2)22C(H(M))C(H(M)),C(h)=(πd+12Γ(h+12)hsin(πh)Γ(2h)Γ(h+d2))12, and the Hurst parameter is a continuous function H:Rd(0,1), where Γ is the usual Gamma function.

  3. The standard fractional Brownian sheets are defined through their covariance function (see ): R(M,M)=12di=1d{|Mi|2Hi+|Mi|2Hi|MiMi|2Hi}, where (H1,,Hd) stands for the multivariate Hurst index in Rd, 0<Hi<1.

  4. The anisotropic fractional Brownian fields are defined through their covariance function (see ): R(M,M)=vH(M)+vH(M)vH(MM), where the variogram vH(x)=22H1γ(H)CH,ϑ1,ϑ2(x)x2H, with H(0,1), γ(H) depends explicitly on H and CH,ϑ1,ϑ2(.) implies incomplete Beta functions and two constants π2ϑ1<ϑ2π2.

The FieldSim package

In the new version 3.2 of the package FieldSim, new features have been added. The most important add is the “process” class and the setProcess function.

An object of class “process” has different slots:

All the examples presented can be defined with the setProcess command (see Table 1). With the following command, the user can set a fBm with Hurst parameter 0.7 on a regular grid of the interval [0,1] (of size 256).

R> linefBm <- setProcess("fBm-line", 0.7)
R> str(linefBm)
Formal class 'process' [package "FieldSim"] with 7 slots
  ..@ name     : chr "fBm"
  ..@ values   : num 0
  ..@ manifold :Formal class 'manifold' [package "FieldSim"] with 4 slots
  .. .. ..@ name    : chr "line"
  .. .. ..@ atlas   : num [1, 1:256] 0 0.00392 0.00784 0.01176 0.01569 ...
  .. .. ..@ distance:function (xi, xj)  
  .. .. ..@ origin  : num [1, 1] 0
  ..@ covf     :function (xi, xj)  
  ..@ parameter: num 0.7
  ..@ values2  : num 0
  ..@ manifold2:Formal class 'manifold' [package "FieldSim"] with 4 slots
  .. .. ..@ name    : chr "line"
  .. .. ..@ atlas   : num [1, 1:256] 0 0.00392 0.00784 0.01176 0.01569 ...
  .. .. ..@ distance:function (xi, xj)  
  .. .. ..@ origin  : num [1, 1] 0

It is worth mentioning that the slot values is empty since there is no simulation done. Then as usual, the user can use the fieldsim function in order to simulate the Gaussian process associated to covf on the manifold grid defined in manifold.

R> fieldsim(linefBm)

In the fieldsim function, we can add the quantity Ne, the number of points of the grid to be simulated in the exact step, and nbNeighbor, the number of neighbors used in the refined step. By default, Ne is equal to the size of the grid given in atlas. The slot values are now set with the simulated values. There exist different visualization procedures to draw the results, for instance:

R> plot(linefbm, "default")

We recall that the discretization grids can be modified with the setAtlas command. Depending on the manifold, there are several types of grids: "regular", "random" and "visualization". For instance,

R> setAtlas(linefBm, "regular", 1000)
R> fieldsim(linefBm)
R> plot(linefBm, "default")

3 The fieldsim procedure for conditional Gaussian fields

In order to build conditional fractional Gaussian fields, we consider a conditioning set N={N1,, Nk}, NiRd, i=1,,k, and the conditioning values x=(x1,,xk)TRk. Then we will say that X~()={X~(M),MRd} is the conditional Gaussian field associated to the field X() (of covariance function R) and to the conditioning pair (N,x) if the finite dimensional laws of X~() is the same as the finite dimensional laws of X() given the event {(X(N1),,X(Nk))T=:XN=x}. We denote by m~() (resp. R~(,)) the mean (resp. covariance) function of the process X~(). The following lemma allows us to determine m~() and R~(,) according to R(,) (sketch of proof is given in ).

Lemma 1. *Let us consider the centered Gaussian vector (Y1,Y2,ZT)TR×R×Rk with the covariance matrix Σ2=(E(Y12)E(Y1Y2)E(Y1ZT)E(Y1Y2)E(Y22)E(Y2ZT)E(ZY1)E(ZY2)E(ZZT)). Suppose that E(ZZT) is invertible. Then the conditional law of (Y1,Y2)T given the event {Z=zRk} is Gaussian with mean

(1)m~=(E(Y1ZT)E(Y2ZT)){E(ZZT)}1z, and covariance matrix (2)Σ~2=(E(Y12)E(Y1Y2)E(Y1Y2)E(Y22))(E(Y1ZT)E(Y2ZT)){E(ZZT)}1(E(ZY1)E(ZY2)).

In the Gaussian field context, Lemma 1 allows us to write down an explicit expression of the mean function and the autocovariance function of the conditional Gaussian field associated to R(,) and to (N,x). Let us put Y1=X(M1) and Y2=X(M2) the values of the field X() at points M1Rd and M2Rd respectively, and Z=XNRk. Therefore, all quantities in (1) and (2) can be expressed in terms of the autocovariance function R. Precisely, E(YiYj)=R(Mi,Mj),(i,j){1,2}2, and E(YiZ)=R(Mi,N),i{1,2},=1,,k. Consequently, the mean function of the conditional Gaussian field is given by

(3)m~(M)=E(X(M)XNT){E(XNXNT)}1x,MRd.

Then the autocovariance function of a conditional Gaussian field (using the (1,2)-coordinate of Equation (2)) is given by (4)R~(M1,M2)=R(M1,M2)E(X(M1)XNT){E(XNXNT)}1E(XNX(M2)). For instance, for k=1, we get m~(M)=R(M,N1)R(N1,N1)x1, and R~(M1,M2)=R(M1,M2)R(M1,N1)R(M2,N1)R(N1,N1).

Let us recall that the goal of this paper is to give a procedure that yields discretization of the sample path of the conditional Gaussian field over a space discretization {Se,Sr} of Rd associated to the n.n.d. autocovariance function R and the conditioning set and values (N,x) . In the sequel, we denote by X~() this sample path. Since the mean function (3) is known, we can consider the centered field X()=X~()m~(). The fieldsim procedure for conditional Gaussian fields proceeds as follows.

Exact simulation step.

Given a space discretization Se, a sample of a centered Gaussian vector (X(M))MSe with covariance matrix R~ given by {R~}i,j=R~(Mi,Mj), Mi,MjSe, is simulated. Here R~ is defined by (4). This simulation is obtained by an algorithm based on Cholesky decomposition of the matrix R~.

Refined simulation step.

Let Sr be the remaining space discretization. For each new point MSr at which we want to simulate the field, X(M) is generated by using only a set of neighbors instead of all the simulated components (as in the accurate simulation step). Precisely, let OM be a neighbors set of M (for the Euclidean distance) and XOM be the space generated by the variables X(M), MOM. Let us remark that the neighbors set is defined with all the already simulated variables (in the accurate and refined simulation step). Let XXOM(M) be the best linear combination of variables of XOM approximating X(M) in the sense that the variance of the innovation εXNM(M)=X(M)XXOM(M), is minimal. The new variable X(M) is obtained by XXOM(M)+Var(εXOM(M))U, where U is a centered and reduced Gaussian variable independent of the already simulated components. Note that the variable XXOM(M) and the variance Var(εXOM(M)) are completely determined by the covariance structure of the sequence X(M), MOM{M}.

Adding the mean.

Finally, we compute X~(M)=X(M)+m~(M) for all M{Se,Sr}.

For storage and computing time, the accurate simulation step must concern only a small number of variables whereas the second step can relate to a larger number of variables. That leads to an effective and fast method to simulate any Gaussian field.

It is worth mentioning that the setProcess command will check if {E(XNXNT)}1 exists for common conditional simulations.

4 Some examples of conditional fractional Gaussian fields

We focus, in this paper, on the conditional Gaussian fields associated to the previously mentioned fields but every other classical Gaussian field can be also simulated: standard bifractional Brownian motion, space-time deformed fractional Brownian motion, etc. (see ). We also consider conditional simulations associated to fractional Gaussian fields on manifolds (hyperboloid and sphere) (see for the covariance function definition).

The procedure fieldsim is extended to the conditional Gaussian fields. We can find the setProcess reference short-card in Table 1.

On the line

The fractional Gaussian processes on the line are fast to simulate.

graphic without alt textgraphic without alt text

Figure 2: Conditional simulations associated to fractional Brownian motion and multifractional Brownian motion. The real time (resp. CPU time) in seconds is equal to 8.430 (resp. 0.043) for the fractional Brownian motion and 14.609 (resp. 0.111) for the multifractional Brownian motion.

Conditional simulations associated to fractional Brownian motion (fBm) and multifractional Brownian motion (mBm) and to the conditioning set N={1/2,3/4,1} and conditioning values x={1,1/2,0} are illustrated on Figure 2. Here the Hurst exponent is H=0.7 for the fBm and H(t)=0.3+0.6t, t[0,1] for the mBm. The processes are simulated on a regular grid of 256 points of [0,1] with only an exact simulation step (Sr=).

They can be obtained with the fieldsim procedure. For instance, the mBm in Figure 2 is obtained with:

R> funcH <- function(x) 0.3 + x * 0.6
R> cond.mBm <- setProcess("cond-mBm-line", 
+    list(Gamma = matrix(c(1/2, 1, 3/4, 0.5, 1, 0), 2, 3), par = funcH))
R> fieldsim(cond.mBm)
R> plot(cond.mBm)

In the simulation below, the points of the set N belong to the visualization grid. When this is not the case, the plot could show a failure for the conditioning in the region of high variability. To avoid this, it is possible to add the points of the set N to the visualization grid. For instance, in the previous example, to add the point 1/6 to the visualization grid, we can use the following lines of code:

R> atlas.cond.mBm <- sort(c(cond.mBm@manifold@atlas[1, ], 1/6))
R> cond.mBm@manifold@atlas <- matrix(atlas.mBm, nrow = 1)

Another solution is to use finer grids which contain the points of the set N.

On the plane

Conditional simulations associated to a fractional Brownian field (for H=0.9) and multifractional Brownian field (for H(t)=0.3+0.6t1) are illustrated in Figure 3. Conditional simulations associated to anisotropic fields (fractional Brownian sheet with H1=0.9, H2=0.3, anisotropic fractional Brownian field with H=0.7, ϑ1=π6 and ϑ2=π3) are presented in Figure 4. For all the fields, we consider the following conditioning set N={(1,k26+1),(k26+1,1),k=0,,26+1}, and conditioning values x=0.

graphic without alt textgraphic without alt text

Figure 3: Conditional simulations associated to a fractional Brownian field (on the left) and a multifractional Brownian field (on the right). The real time (resp. CPU time) in seconds is equal to 1270.911 (resp. 6.769) for the fractional Brownian field and 1782.533 (resp. 9.953) for the multifractional Brownian field.

graphic without alt textgraphic without alt text

Figure 4: Conditional simulations associated to a fractional Brownian sheet (on the left) and an anisotropic fractional Brownian field. The real time (resp. CPU time) in seconds is equal to 728.605 (resp. 4.223) for the fractional Brownian sheet and 2995.644 (resp. 14.782) for the anisotropic fractional Brownian field.

All the processes are simulated on a regular grid of 4096 points of [0,1]2 with 100 points for the exact simulation step and 3996 for the refined step (with 4 neighbors). For instance, the conditional Gaussian field associated to anistropic fractional Brownian field on [0,1]2 (see Figure 4) is given by

R> Ng <- 2^6 + 1
R> x <- seq(from = 0, to = 1, length = Ng)
R> G <- cbind(rbind(rep(1, Ng - 1), x[2:Ng], rep(0, Ng - 1)), 
+    rbind(x[2:(Ng - 1)], rep(1, Ng - 2), rep(0, Ng - 2)))  
R> condfBm2d <- setProcess("cond-afBf-plane", 
+    list(Gamma = G, par = list(H = 0.7, theta1 = pi/6, theta2 = pi/3)))
R> setAtlas(condfBm2d, "visualization", 6)
R> fieldsim(condfBm2d, Ne = 100, nbNeighbor = 4)
R> plot(condfBm2d, theta = 120, phi = 30, expand = 0.5)

It is worth emphasizing that, for a fixed size of the simulation grid, the simulation time of the fieldsim procedure depends on the number of conditioning points (see for instance the fractional Brownian field on Figures 3 and 6 for 129 and 39 conditioning points respectively). But this variation is small compared to the variation due to the size of the simulation grid.

On the hyperboloid and on the sphere

Conditional simulations can be extended to fractional Gaussian fields on manifolds associated to the fractional Brownian field on the hyperboloid with H=0.7, N={(0,1,2),(0,2,5)},x=(5,5), and a conditional fractional Brownian field on the sphere with H=0.4, N={(0,0,1),(12,0,32)},x=(5,5). The two processes are simulated on a regular grid of 5400 points of R3 with 100 points for the exact simulation step and 5300 for the refined step (with 4 neighbors).

graphic without alt textgraphic without alt text

Figure 5: Conditional simulations associated to the fractional Brownian field on the hyperboloid and on the sphere. The real time (resp. CPU time) in seconds is equal to 54.567 (resp. 0.293) for the hyperboloid and 188.807 (resp. 14.216) for the sphere.

The conditional simulations associated to the fractional Brownian field on the sphere (see Figure 5) are obtained with

R> Gamma <- matrix(c(0, 0, 1, 5, 0.5, 0, sqrt(3)/2, -5), 4, 2)
R> sphere.cond.fBm <- setProcess("cond-fBm-sphere", list(Gamma = Gamma, par = 0.4))
R> setAtlas(sphere.cond.fBm, "visualization", 30) 
R> fieldsim(sphere.cond.fBm, Ne = 100, nbNeighbor = 4) 
R> plot(sphere.cond.fBm)

5 Conclusion and perspectives

We propose a generic method to simulate multidimensional conditional fractional Gaussian fields.

Our method is valid for any Gaussian field and on any (non regular) grid of points as soon as the covariance function is available. This method is constructed to be universal (conditional simulation, simulation on a manifold) and is consecutively not as fast as other methods defined for specific fields. In the near future, the FieldSim package should also possess such specific methods.

Our method is adapted to conditional simulations and, consequently, permits now to simulate easily several natural scenes (clouds, mountains) with valleys and fixed topographic points. Such a simulation is presented in Figure 6.

graphic without alt text
Figure 6: Natural scene simulation. Here a 65×65 regular grid fractional Brownian field of Hurst parameter H=0.8 with 39 conditioning points (in red) is simulated. The real time (resp. CPU time) in seconds is equal to 563.754 (resp. 3.550).
Table 1: The “process” class. "fBm" for fractional Brownian motion, "mBm" for multifractional Brownian motion, "2pfBm" for the standard bi-fractional Brownian motion, "stdfBm" for the space-time deformed fractional Brownian motion, "AfBf" for anisotropic fractional Brownian field and "fBs" for fractional Brownian sheet.
name of the process parameter
On [0,1]
"fBm-line" numeric
"mBm-line" function
"2pfBm-line" list(H = numeric, K = numeric)
"stdfBm-line" list(H = numeric, sigma = function, tau = function)
"cond-fBm-line" list(Gamma = matrix, par = numeric)
"cond-mBm-line" list(Gamma = matrix, par = function)
"cond-2pfBm-line" list(Gamma = matrix, par = list(H = numeric, K = numeric))
On the square [0,1]2
"fBm-plane" numeric
"mBm-plane" function
"2pfBm-plane" list(H = numeric, K = numeric)
"afBf-plane" list(H = numeric, theta1 = numeric, theta2 = numeric)
"fBs-plane" vector
"cond-fBm-plane" list(Gamma = matrix, par = numeric)
"cond-mBm-line" list(Gamma = matrix, par = function)
"cond-afBf-plane" list(Gamma = matrix, par = list(H = numeric, theta1 = numeric,
theta2 = numeric))
"cond-fBs-line" list(Gamma = matrix, par = vector)
On the hyperboloid
"fBm-hyperboloid" numeric
"cond-fBm-hyperboloid" list(Gamma = matrix, par = numeric)

CRAN packages used

FieldSim, RandomFields

CRAN Task Views implied by cited packages

Note

This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.

Footnotes

    References

    A. Benassi, S. Jaffard and D. Roux. Elliptic Gaussian random processes. Revista Mathematica Iberoamericana, 18: 19–89, 1997.
    A. Bonami and A. Estrade. Anisotropic analysis of some Gaussian models. Journal of Fourier Analysis and Applications, 9(3): 215–236, 2003.
    A. Brouste, J. Istas and S. Lambert-Lacroix. On fractional Gaussian random fields simulations. Journal of Statistical Software, 23(1): 1–23, 2007. DOI 10.18637/jss.v023.i01.
    A. Brouste, J. Istas and S. Lambert-Lacroix. On simulation of manifold indexed fractional Gaussian fields. Journal of Statistical Software, 34(4): 1–14, 2010. DOI 10.18637/jss.v034.i04.
    A. Brouste and S. Lambert-Lacroix. FieldSim: Random fields (and bridges) simulations. 2015. URL https://CRAN.R-project.org/package=FieldSim. R package version 3.2.1.
    J. F. Coeurjolly. Simulation and identification of the fractional Brownian motion: A bibliographical and comparative study. Journal of Statistical Software, 5(7): 1–53, 2000. DOI 10.18637/jss.v005.i07.
    Z. Gelbaum and M. Titus. Simulation of fractional Brownian surfaces via spectral synthesis on manifolds. IEEE Transactions on Image Processing, 23(10): 4383–4388, 2014.
    A. Kamont. On the fractional anisotropic Wiener fields. Journal of Probability and Mathematical Statistics, 18: 85–98, 1996.
    A. Kolmogorov. Wienersche Spiralen und einige andere interessante Kurven im Hilbertschen Raum. (German). C. R. (Doklady) Academy of Sciences URSS, 26: 115–118, 1940.
    B. Mandelbrot and J. Van Ness. Fractional Brownian motions, fractional noises and application. SIAM Review, 10: 422–437, 1968.
    H. O. Peitgen and D. Saupe. The science of fractal images. Springer-Verlag, 1988.
    R. Peltier and J. Levy-Véhel. Multifractional Brownian motion: Definition and preliminary results. RR 2645. INRIA. 1996. URL http://hal.inria.fr/inria-00074045/fr/.
    V. Piterbag. Asymptotic methods in the theory of gaussian processes and fields. American Mathematical Society, 1996.
    G. Samorodnitsky and M. Taqqu. Stable non-Gaussian random processes: Stochastic models with infinite variance. New York: Chapman & Hall, 1994.
    M. Schlather, A. Malinowski, P. J. Menck, M. Oesting and K. Strokorb. Analysis, simulation and prediction of multivariate random fields with package RandomFields. Journal of Statistical Software, 63(8): 1–25, 2015. DOI 10.18637/jss.v063.i08.
    M. Schlather, A. Malinowski, M. Oesting, D. Boecker, K. Strokorb, S. Engelke, J. Martini, F. Ballani, O. Moreva, P. J. Menck, et al. RandomFields: Simulation and analysis of random fields. 2016. URL http://CRAN.R-project.org/package=RandomFields. R package version 3.1.8.
    M. Stein. Fast and exact simulation of fractional Brownian surfaces. Journal of Computational and Graphical Statistics, 11(3): 587–599, 2002.
    A. Wood and G. Chan. Simulation of stationary Gaussian processes in [0,1]d. Journal of Computational and Graphical Statistics, 3(4): 409–432, 1994.
    Z. Yin. New methods for simulation of fraction Brownian motion. Journal of Computational Physics, 127: 66–72, 1996.

    Reuse

    Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".

    Citation

    For attribution, please cite this work as

    Brouste, et al., "Conditional Fractional Gaussian Fields with the Package FieldSim", The R Journal, 2016

    BibTeX citation

    @article{RJ-2016-003,
      author = {Brouste, Alexandre and Istas, Jacques and Lambert-Lacroix, Sophie},
      title = {Conditional Fractional Gaussian Fields with the Package FieldSim},
      journal = {The R Journal},
      year = {2016},
      note = {https://rjournal.github.io/},
      volume = {8},
      issue = {1},
      issn = {2073-4859},
      pages = {38-47}
    }