RandomFields (version 3.1.12)

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, distances, dim, likelihood,
                estimate_variance =NA, ...)

Arguments

model

object of class RMmodel; the covariance or variogram model, which is to be evaluated

x

vector or $(n \times \code{dim})$-matrix, where $n$ is the number of points at which the covariance function is to be evaluated; in particular, if the model is isotropic or dim=1 then x is a vector. x

y

second vector or matrix for non-stationary covariance functions

z

z-component of point if xyzT-specification of points is used

T

T-component of point if xyzT-specification of points is used

grid

boolean; whether xyzT specify a grid

data

vector or matrix of values measured at coord; If a matrix is given then the columns are interpreted as independent realisations. If also a time component is given, then in the data the indices for the spatial components run the fastest.

distances

vector; the lower triangular part of the distance matrix column-wise; equivalently the upper triangular part of the distance matrix row-wise; either x or distances must be missing

dim

dimension of the coordinate space in which the model is applied; only necesary for given distances

likelihood

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

estimate_variance

logical or NA. See Details.

...

for advanced further options and control arguments for the simulation that are passed to and processed by RFoptions

Value

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

Details

The function calculates the likelihood for data of a Gassian 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.

See Also

Bayesian, RMmodel, RFfit, RFsimulate, RFlinearpart.

Examples

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


require("mvtnorm")

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

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

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

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

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

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


pts <- c(4, 5)
repet <- c(2, 3)
trend <- 2 * sin(R.p(new="isotropic")) + 3
model <- 2 * RMexp() + trend
x <- y <- data <- 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))
  data[[i]] <- as.matrix(RFsimulate(model, x=x[[i]]$x, y=x[[i]]$y,
                                    n=repet[i], spC = FALSE))
}

print(unix.time(likeli <- RFlikelihood(model, x, data=data)))
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(data[[p]])) {
    print(unix.time(dn <- dmvnorm(data[[p]][,i], mean=tr, sigma=C,
log=TRUE)))
    L <- L + dn
  }
}
print(L)


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



FinalizeExample()