1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
|
\name{RFloglikelihood}
\alias{RFloglikelihood}
\alias{RFlikelihood}
\title{Likelihood and estimation of linear models}
\description{
\command{\link{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{
\item{model,params}{\argModel}
\item{x}{\argX}
\item{y,z}{\argYz}
\item{T}{\argT}
\item{grid}{\argGrid}
\item{distances,dim}{\argDistances}
\item{data}{\argData}
%\item{set}{integer. See section Value for details.}
\item{likelihood}{ Not programmed yet. Character.
Choice of kind of likelihood ("full", "composite", etc.),
see also \code{likelihood} for \command{\link{RFfit}}
in \command{\link{RFoptions}}.
}
%\item{log}{logical. If \code{TRUE} the loglikelihood is returned.
% }
\item{estimate_variance}{logical or \code{NA}. See Details.
}
\item{...}{\argDots}
}
\details{
The function calculates the likelihood for data of a Gaussian process
with given covariance structure.
The covariance structure may not have \code{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{
\command{\link{RFloglikelihood}} returns a list
containing the likelihood, the log likelihood, and
the global variance (if estimated -- see details).
}
\me
\seealso{
\link{Bayesian},
\command{\link{RMmodel}},
\command{\link{RFfit}},
\command{\link{RFsimulate}},
\command{\link{RFlinearpart}}.
}
\examples{\dontshow{StartExample()}
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))
\dontshow{FinalizeExample()}}
\keyword{spatial}
|
