RandomFields: RFloglikelihood – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

RFloglikelihood

Likelihood and estimation of linear models

Description

RFloglikelihood returns the log likelihood for Gaussian random fields. In case NAs are given that refer to linear modeling, the ML of the linear model is returned.

Usage

RFlikelihood(model, x, y = NULL, z = NULL, T = NULL, grid = NULL,
                data, params, distances, dim, likelihood,
                estimate_variance =NA, ...)

Arguments

`model,params`	object of class `RMmodel`, `RFformula` or `formula`; best is to consider the examples below, first. The argument `params` is a list that specifies free parameters in a formula description, see RMformula.
`x`	vector of x coordinates, or object of class `GridTopology` or `raster`; for more options see RFsimulateAdvanced.
`y,z`	optional vectors of y (z) coordinates, which should not be given if `x` is a matrix.
`T`	optional vector of time coordinates, `T` must always be an equidistant vector. Instead of `T=seq(from=From, by=By, len=Len)`, one may also write `T=c(From, By, Len)`.
`grid`	logical; the function finds itself the correct value in nearly all cases, so that usually `grid` need not be given. See also RFsimulateAdvanced.
`distances,dim`	another alternative for the argument `x` to pass the (relative) coordinates, see RFsimulateAdvanced.
`data`	matrix, data.frame or object of class `RFsp`; If a matrix is given the ordering of the colums is the following: space, time, multivariate, repetitions, i.e. the index for the space runs the fastest and that for repetitions the slowest.

likelihood

Not programmed yet. Character. Choice of kind of likelihood ("full", "composite", etc.), see also likelihood for RFfit in RFoptions.

`estimate_variance`	logical or `NA`. See Details.
`...`	for advanced use: further options and control arguments for the simulation that are passed to and processed by `RFoptions`. If `params` is given, then `...` may include also the variables used in `params`.

Details

The function calculates the likelihood for data of a Gaussian process with given covariance structure. The covariance structure may not have NA values in the parameters except for a global variance. In this case the variance is returned that maximizes the likelihood. Additional to the covariance structure the model may include a trend. The latter may contain unknown linear parameters. In this case again, the unknown parameters are estimated, and returned.

Value

RFloglikelihood returns a list containing the likelihood, the log likelihood, and the global variance (if estimated – see details).

Author(s)

Martin Schlather, schlather@math.uni-mannheim.de, https://www.wim.uni-mannheim.de/schlather/

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again


requireNamespace("mvtnorm")

pts <- 4
repet <- 3
model <- RMexp()
x <- runif(n=pts, min=-1, max=1)
y <- runif(n=pts, min=-1, max=1)
dta <- as.matrix(RFsimulate(model, x=x, y=y, n=repet, spC = FALSE))
print(cbind(x, y, dta))
print(system.time(likeli <- RFlikelihood(model, x, y, data=dta)))
str(likeli, digits=8)

L <- 0
C <- RFcovmatrix(model, x, y)
for (i in 1:ncol(dta)) {
  print(system.time(dn <- mvtnorm::dmvnorm(dta[,i], mean=rep(0, nrow(dta)),
sigma=C, log=TRUE)))
  L <- L + dn
}
print(L)
stopifnot(all.equal(likeli$log, L))






pts <- 4
repet <- 1
trend <- 2 * sin(R.p(new="isotropic")) + 3
#trend <- RMtrend(mean=0)
model <- 2 * RMexp() + trend
x <- seq(0, pi, len=pts)
dta <- as.matrix(RFsimulate(model, x=x, n=repet, spC = FALSE))
print(cbind(x, dta))

print(system.time(likeli <- RFlikelihood(model, x, data=dta)))
str(likeli, digits=8)

L <- 0
tr <- RFfctn(trend, x=x, spC = FALSE)
C <- RFcovmatrix(model, x)
for (i in 1:ncol(dta)) {
  print(system.time(dn <- mvtnorm::dmvnorm(dta[,i], mean=tr, sigma=C,log=TRUE)))
  L <- L + dn
}
print(L)

stopifnot(all.equal(likeli$log, L))







pts <- c(3, 4)
repet <- c(2, 3)
trend <- 2 * sin(R.p(new="isotropic")) + 3
model <- 2 * RMexp() + trend
x <- y <- dta <- list()
for (i in 1:length(pts)) {
  x[[i]] <- list(x = runif(n=pts[i], min=-1, max=1),
                 y = runif(n=pts[i], min=-1, max=1))
  dta[[i]] <- as.matrix(RFsimulate(model, x=x[[i]]$x, y=x[[i]]$y,
                                    n=repet[i], spC = FALSE))
}

print(system.time(likeli <- RFlikelihood(model, x, data=dta)))
str(likeli, digits=8)

L <- 0
for (p in 1:length(pts)) {
  tr <- RFfctn(trend, x=x[[p]]$x, y=x[[p]]$y,spC = FALSE)
  C <- RFcovmatrix(model, x=x[[p]]$x, y=x[[p]]$y)
  for (i in 1:ncol(dta[[p]])) {
    print(system.time(dn <- mvtnorm::dmvnorm(dta[[p]][,i], mean=tr, sigma=C,
                                    log=TRUE)))
    L <- L + dn
  }
}
print(L)
stopifnot(all.equal(likeli$log, L))

RandomFields

Simulation and Analysis of Random Fields

v3.3.10

GPL (>= 3)

Authors

Martin Schlather [aut, cre], Alexander Malinowski [aut], Marco Oesting [aut], Daphne Boecker [aut], Kirstin Strokorb [aut], Sebastian Engelke [aut], Johannes Martini [aut], Felix Ballani [aut], Olga Moreva [aut], Jonas Auel[ctr], Peter Menck [ctr], Sebastian Gross [ctr], Ulrike Ober [ctb], Paulo Ribeiro [ctb], Brian D. Ripley [ctb], Richard Singleton [ctb], Ben Pfaff [ctb], R Core Team [ctb]

Initial release