Learn R Programming

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

Run this code
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()

Run the code above in your browser using DataLab