This notebook presents sevral possibilities of the SPDE procedure
Some global constants
ndat = 1000
rangev = 0.2
sill = 1.
nugget = 0.1
law_set_random_seed(123)
## NULL
Creating the Data Base for conditioning information
dat = Db_create()
dat["x"] = VectorHelper_simulateUniform(ndat)
dat["y"] = VectorHelper_simulateUniform(ndat)
dat$setLocators(c("x","y"),ELoc_X())
## NULL
dat
##
## Data Base Characteristics
## =========================
##
## Data Base Summary
## -----------------
## File is organized as a set of isolated points
## Space dimension = 2
## Number of Columns = 2
## Maximum Number of UIDs = 2
## Total number of samples = 1000
##
## Variables
## ---------
## Column = 0 - Name = x - Locator = x1
## Column = 1 - Name = y - Locator = x2
Creating the Data Base for the output grid
grid = DbGrid_create(nx=c(50,50),dx=c(0.02,0.02))
Creating the Meshing (Turbo) based on extended grid
gridExt = DbGrid_create(nx=c(75,75),dx=c(0.02,0.02),x0=c(-0.25,-0.25))
mesh = MeshETurbo(gridExt)
mesh
##
## Grid characteristics:
## ---------------------
## Origin : -0.250 -0.250
## Mesh : 0.020 0.020
## Number : 75 75
##
## Turbo Meshing
## =============
## Euclidean Geometry
## Space Dimension = 2
## Number of Apices per Mesh = 3
## Number of Meshes = 10952
## Number of Apices = 5625
##
## Bounding Box Extension
## ----------------------
## Dim #1 - Min:-0.25 - Max:1.23
## Dim #2 - Min:-0.25 - Max:1.23
Creating the Model
model = Model_createFromParam(type=ECov_BESSEL_K(),param=1,range=rangev,sill=sill)
model$addCovFromParam(ECov_NUGGET(),sill=nugget)
## NULL
model
##
## Model characteristics
## =====================
## Space dimension = 2
## Number of variable(s) = 1
## Number of basic structure(s) = 2
## Number of drift function(s) = 0
## Number of drift equation(s) = 0
##
## Covariance Part
## ---------------
## K-Bessel (Third Parameter = 1)
## - Sill = 1.000
## - Range = 0.200
## - Theo. Range = 0.058
## Nugget Effect
## - Sill = 0.100
## Total Sill = 1.100
Building the SPDE environment for Non-conditional environments
spdeS = SPDE()
spdeS$init(model=model, field=grid, calc=ESPDECalcMode_SIMUNONCOND(), mesh=mesh)
## NULL
spdeS$compute()
## NULL
Apply the SPDE system in order to obtain a non-conditional simulation on the input and output Data Bases
iuid = spdeS$query(dat)
iuid = spdeS$query(grid)
Representing the resulting non-conditional simulation on the Grid
p = ggplot()
p = p + plot.grid(grid, "spde.simu")
p = p + plot.decoration(title="Non Conditional Simulation")
ggPrint(p)
Preparing the SPDE environment
spdeK = SPDE()
spdeK$init(model=model, field=grid, data=dat, calc=ESPDECalcMode_KRIGING(), mesh=mesh)
## NULL
spdeK$compute()
## NULL
Apply the SPDE system in order to obtain an estimation on the output Data Base
iuid = spdeK$query(grid, namconv=NamingConvention("spde",FALSE))
Representing the resulting estimation on the Grid
p = ggplot()
p = p + plot.grid(grid, "spde.kriging")
p = p + plot.decoration(title="Estimation")
ggPrint(p)
Extracting the internal elements from 'spdeK' in order to perform calculations by hand (if necessary).
Q = PrecisionOpCs(mesh, model)$getQToTriplet(flag_from_1=TRUE)$toTL()
Aproj = ProjMatrix(dat,mesh)$getAprojToTriplet(flag_from_1=TRUE)$toTL()
Posterior calculations
size = dim(Q)[1]
cholQ = Cholesky(Q)
u = VectorHelper_simulateGaussian(size)
Qop = Q + 1/nugget * t(Aproj) %*% Aproj