RMtrend: Trend Model

Description

RMtrend is a pure trend model with covariance 0.

Usage

RMtrend(mean, plane, polydeg, polycoeff, arbitraryfct, fctcoeff)

Arguments

mean

optional; should be a vector of length $p$, where $p$ is the number of variables taken into account by the corresponding multivariate random field $(Z_1(\cdot),\ldots,Z_p(\cdot))$; the $i$-th component of mean is interpreted as constant

plane

optional; should be a $d \times p$-matrix where $d$ is the dimension of the random field and $p$ the number of variables considered; corresponds to a trend described by $p$ hyperplanes. The mean of $Z_i(x)$ with $x= (x_1,\ldots,x_d)$ is given by $$pl

polydeg

optional; should be an integer vector of length $p$; indicates that the mean of $Z_i(\cdot)$ is given by a multivariate polynomial of degree less than or equal to the $i$-th component of polydeg; the coefficients are either assumed to be

polycoeff

optional; should be a vector of length

(\binom{p o l y d e g_{1} + d}{d}) + \dots + (\binom{p o l y d e g_{p} + d}{d})

which is the number of monomial basis functions up to degree polydeg_1, ..., polydeg_p in a $d$-dimensional space. Each compo

arbitraryfct

optional; should be a $p$-variate function; the $i$-th component of arbitraryfct describes the trend surface of $Z_i(\cdot)$; the arguments of this function should be location (and time) corresponding to the random field to be modelled. I

fctcoeff

optional; should be numerical; determines the coefficient belonging to arbitraryfct, i.e. the trend is given by fctcoeff * arbitraryfct. Note that the coefficient is the same for each component of arbitraryfct;

Value

RMtrend returns an object of class RMmodel.

Details

If different trend arguments are given, the corresponding trend components are added. Equivalently, + (see RFformula) or RMplus can be used to add trend terms.

Note that this function refers to trend surfaces in the geostatistical framework. Fixed effects in the mixed models framework are implemented, as well (see RFformula).

The order of the monomial basis functions and the corresponding coefficients given by polycoeff is the following: The first ${polydeg_1+d\choose d}$ components belong to the trend polynomial of the first variable, the following ${polydeg_2+d\choose d}$ ones to the second one, and so on. Within one trend polynomial the monomial basis functions are ordered by the powers in an ascending way such that the power of the first component varies fastest; e.g. the monomial basis functions up to degree $k$ in a two-dimensional space are given by $1, x, \dots, x^{k}, y, x y, \dots, x^{k - 1} y, y^{2}, \dots, x y^{k - 1}, y^{k} .$

References

Chiles, J. P., Delfiner, P. (1999) Geostatistics: Modelling Spatial Uncertainty. New York: John Wiley & Sons.

Examples

Run this code

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again
##################################################
# Example 1: # 
# Simulate from model with a plane trend surface #
##################################################

#trend: 1 + x - y, cov: exponential with variance 0.01
model <- ~ RMtrend(mean=1, plane = c(1,-1)) + RMexp(var=0.01)
#equivalent model:
model <- ~ RMtrend(polydeg=1,polycoeff=c(1,1,-1)) + RMexp(var=0.01)
#Simulation
x <- 0:10
simulated <- RFsimulate(model=model, x=x, y=x, grid=TRUE)
plot(simulated)


####################################################################
#
# Example 2: Simulate and fit a multivariate geostatistical model
#
####################################################################
 
# Simulate a bivariate Gaussian random field with trend
# m_1((x,y)) = x + 2*y and m_2((x,y)) = 3*x + 4*y
# where m_1 is a hyperplane describing the trend for the first response
# variable and m_2 is the trend for the second one;
# the covariance function is the sum of a bivariate Whittle-Matern model
# and a multivariate nugget effect
x <- y <- 0:10
model <- RMtrend(plane=matrix(c(1,2,3,4), ncol=2)) + 
         RMparswm(nu=c(1,1)) + RMnugget(var=0.5)
multi.simulated <- RFsimulate(model=model, x=x, y=y, grid=TRUE)
plot(multi.simulated)

\dontrun{
# Fit the Gaussian random field with unknown trend for the second
# response variable and unknown variances
model.na <- RMtrend(plane=matrix(c(1, 2, NA, NA), ncol=2)) + 
            RMparswm(nu=c(1,1), var=NA) + RMnugget(var=NA)
fit <- RFfit(model=model.na, data=multi.simulated)
}

\dontrun{
##################################################
#
# Example 3:  Simulation and estimation for model with 
#             arbitrary trend functions
#
##################################################

#Simulation
# trend: 2*sin(x) + 0.5*cos(y), cov: spherical with scale 3
model <- ~ RMtrend(arbitraryfct=function(x) sin(x),
 fctcoeff=2) +
 RMtrend(arbitraryfct=function(y) cos(y),
 fctcoeff=0.5) +
 RMspheric(scale=3)
x <- seq(-4*pi, 4*pi, pi/10)
simulated <- RFsimulate(model=model, x=x, y=x, grid=TRUE)
plot(simulated)

################# ?? !!
#Estimation, part 1
# estimate coefficients and scale
model.est <- ~ RMtrend(arbitraryfct=function(x) sin(x), fctcoeff=1) +
 RMtrend(arbitraryfct=function(y) cos(y), fctcoeff=1) +
 RMspheric(scale=NA)
estimated <- RFfit(model=model.est, x=x, y=x, grid=TRUE,
 data=simulated@data, mle.methods="ml")


#Estimation
# estimate coefficients and scale
model.est <- ~ RMtrend(arbitraryfct=function(x) sin(x)) +
 RMtrend(arbitraryfct=function(y) cos(y)) +
 RMspheric(scale=NA)
estimated <- RFfit(model=model.est, x=x, y=x, grid=TRUE,
 data=simulated@data, mle.methods="ml")



##################################################
#
# Example 4: Simulation and estimation for model with 
#            polynomial trend 
#
##################################################

#Simulation
# trend: 2*x^2 - 3 y^2, cov: whittle-matern with nu=1,scale=0.5
model <- ~ RMtrend(arbitraryfct=function(x) 2*x^2 - 3*y^2,
 fctcoeff=1) + RMwhittle(nu=1, scale=0.5)
# equivalent model:
model <- ~ RMtrend(polydeg=2, polycoeff=c(0,0,2,0,0,-3))
x <- 0:20		 
simulated <- RFsimulate(model=model, x=x, y=x, grid=TRUE)
plot(simulated)

#Estimation
# estimate nu and the trend term assuming that it is a polynomial
# of degree 2
model.est <- ~ RMtrend(polydeg=2) + RMwhittle(nu=NA, scale=0.5)
estimated <- RFfit(model=model.est, x=x, y=x, grid=TRUE,
 data=simulated@data, mle.methods="ml")
}
FinalizeExample()

Run the code above in your browser using DataLab

Data Engineering and BI courses are free this week!