# R/Matrices.R In LaplacesDemon: Complete Environment for Bayesian Inference

```###########################################################################
# Matrices                                                                #
#                                                                         #
# These are utility functions for matrices.                               #
###########################################################################

as.indicator.matrix <- function(x)
{
n <- length(x)
x <- as.factor(x)
X <- matrix(0, n, length(levels(x)))
X[(1:n) + n*(unclass(x)-1)] <- 1
dimnames(X) <- list(names(x), levels(x))
return(X)
}
as.inverse <- function(x)
{
if(!is.matrix(x)) x <- matrix(x)
if(!is.square.matrix(x)) stop("x must be a square matrix.")
if(!is.symmetric.matrix(x)) stop("x must be a symmetric matrix.")
tol <- .Machine\$double.eps
options(show.error.messages=FALSE)
xinv <- try(solve(x))
if(inherits(xinv, "try-error")) {
k <- nrow(x)
eigs <- eigen(x, symmetric=TRUE)
if(min(eigs\$values) < tol) {
tolmat <- diag(k)
for (i in 1:k)
if(eigs\$values[i] < tol) tolmat[i,i] <- 1/tol
else tolmat[i,i] <- 1/eigs\$values[i]
}
else tolmat <- diag(1/eigs\$values, nrow=length(eigs\$values))
xinv <- eigs\$vectors %*% tolmat %*% t(eigs\$vectors)
}
options(show.error.messages=TRUE)
xinv <- as.symmetric.matrix(xinv)
return(xinv)
}
as.parm.matrix <- function(x, k, parm, Data, a=-Inf, b=Inf, restrict=FALSE,
chol=FALSE)
{
X <- matrix(0, k, k)
if(restrict == TRUE) {
X[upper.tri(X, diag=TRUE)] <- c(1,
parm[grep(deparse(substitute(x)),
Data[["parm.names"]])])}
else {
X[upper.tri(X, diag=TRUE)] <- parm[grep(deparse(substitute(x)),
Data[["parm.names"]])]}
if(chol == TRUE) {
if(a != -Inf | b != Inf) {
x <- as.vector(X[upper.tri(X, diag=TRUE)])
x.num <- which(x < a)
x[x.num] <- a
x.num <- which(x > b)
x[x.num] <- b
X[upper.tri(X, diag=TRUE)] <- x
diag(X) <- abs(diag(X))
}
X[lower.tri(X)] <- 0
return(X)
}
X[lower.tri(X)] <- t(X)[lower.tri(X)]
if(a != -Inf | b != Inf) {
x <- as.vector(X[upper.tri(X, diag=TRUE)])
x.num <- which(x < a)
x[x.num] <- a
x.num <- which(x > b)
x[x.num] <- b
X[upper.tri(X, diag=TRUE)] <- x
X[lower.tri(X)] <- t(X)[lower.tri(X)]
}
if(!is.symmetric.matrix(X)) X <- as.symmetric.matrix(X)
if(!exists("LDEnv")) LDEnv <- new.env()
if(restrict == FALSE) {
if(is.positive.definite(X)) {
assign("LaplacesDemonMatrix", as.vector(X[upper.tri(X,
diag=TRUE)]), envir=LDEnv)}
else {
if(exists("LaplacesDemonMatrix", envir=LDEnv)) {
X[upper.tri(X,
diag=TRUE)] <- as.vector(get("LaplacesDemonMatrix",
envir=LDEnv))
X[lower.tri(X)] <- t(X)[lower.tri(X)]}
else {X <- diag(k)}}
}
if(restrict == TRUE) {
if(is.positive.definite(X)) {
assign("LaplacesDemonMatrix", as.vector(X[upper.tri(X,
diag=TRUE)][-1]), envir=LDEnv)}
else {
if(exists("LaplacesDemonMatrix", envir=LDEnv)) {
X[upper.tri(X, diag=TRUE)] <- c(1,
as.vector(get("LaplacesDemonMatrix",
envir=LDEnv)))
X[lower.tri(X)] <- t(X)[lower.tri(X)]
if(!is.symmetric.matrix(X)) X <- as.symmetric.matrix(X)
}
else {X <- diag(k)}}
}
return(X)
}
as.positive.definite <- function(x)
{
eig.tol <- 1e-06
conv.tol <- 1e-07
posd.tol <- 1e-08
iter <- 0; maxit <- 100
n <- ncol(x)
D_S <- x
D_S[] <- 0
X <- x
converged <- FALSE
conv <- Inf
while (iter < maxit && !converged) {
Y <- X
R <- Y - D_S
e <- eigen(R, symmetric=TRUE)
Q <- e\$vectors
d <- e\$values
p <- d > eig.tol * d[1]
if(!any(p))
stop("Matrix seems negative semi-definite.")
Q <- Q[, p, drop=FALSE]
X <- tcrossprod(Q * rep(d[p], each=nrow(Q)), Q)
D_S <- X - R
conv <- norm(Y - X, "I") / norm(Y, "I")
iter <- iter + 1
converged <- (conv <= conv.tol)
}
if(!converged) {
warning("as.positive.definite did not converge in ", iter,
" iterations.")}
e <- eigen(X, symmetric=TRUE)
d <- e\$values
Eps <- posd.tol * abs(d[1])
if(d[n] < Eps) {
d[d < Eps] <- Eps
Q <- e\$vectors
o.diag <- diag(X)
X <- Q %*% (d * t(Q))
D <- sqrt(pmax(Eps, o.diag)/diag(X))
X[] <- D * X * rep(D, each=n)}
X <- as.symmetric.matrix(X)
return(X)
}
as.positive.semidefinite <- function(x)
{
if(!is.matrix(x)) x <- matrix(x)
if(!is.square.matrix(x)) stop("x must be a square matrix.")
if(!is.symmetric.matrix(x)) stop("x must be a symmetric matrix.")
iter <- 0; maxit <- 100
converged <- FALSE
while (iter < maxit && !converged) {
iter <- iter + 1
out <- eigen(x=x, symmetric=TRUE)
mGamma <- t(out\$vectors)
vLambda <- out\$values
vLambda[vLambda < 0] <- 0
x <- t(mGamma) %*% diag(vLambda) %*% mGamma
x <- as.symmetric.matrix(x)
if(is.positive.semidefinite(x)) converged <- TRUE
}
if(converged == FALSE) {
warning("as.positive.semidefinite did not converge in ", iter,
" iterations.")}
return(x)
}
as.symmetric.matrix <- function(x, k=NULL)
{
if(is.vector(x)) {
if(any(!is.finite(x))) stop("x must have finite values.")
if(is.null(k)) k <- (-1 + sqrt(1 + 8 * length(x))) / 2
symm <- matrix(0, k, k)
symm[lower.tri(symm, diag=TRUE)] <- x
symm2 <- symm
symm2[upper.tri(symm2, diag=TRUE)] <- 0
symm <- symm + t(symm2)
}
else if(is.matrix(x)) {
if(!is.square.matrix(x)) stop("x must be a square matrix.")
if(any(!is.finite(diag(x))))
stop("The diagonal of x must have finite values.")
symm <- x
x.lower.fin <- FALSE; x.upper.fin <- FALSE
if(all(is.finite(x[lower.tri(x, diag=TRUE)])))
x.lower.fin <- TRUE
if(all(is.finite(x[upper.tri(x, diag=TRUE)])))
x.upper.fin <- TRUE
if(x.lower.fin) symm[upper.tri(x)] <- t(x)[upper.tri(x)]
else if(x.upper.fin) symm[lower.tri(x)] <- t(x)[lower.tri(x)]
else {
new.up <- x[upper.tri(x)]
new.low <- x[lower.tri(x)]
new.up[which(!is.finite(new.up))] <- t(x)[lower.tri(x)][which(!is.finite(new.up))]
new.low[which(!is.finite(new.low))] <- t(x)[upper.tri(x)][which(!is.finite(new.low))]
if(any(!is.finite(c(new.up, new.low))))
stop("Off-diagonals in x must have finite values.")
else {
symm[upper.tri(symm)] <- new.up
symm[lower.tri(symm)] <- new.low
}
}
}
else stop("x must be a vector or matrix.")
return(symm)
}
.colVars <- function(X)
{
N <- nrow(X)
Y <- X - matrix(colMeans(X), N, ncol(X), byrow=TRUE)
Z <- colMeans(Y*Y)*N/{N-1}
return(Z)
}
Cov2Cor <- function(Sigma)
{
if(missing(Sigma)) stop("Sigma is a required argument.")
if(any(!is.finite(Sigma))) stop("Sigma must have finite values.")
if(is.matrix(Sigma)) {
if(!is.positive.definite(Sigma))
stop("Sigma is not positive-definite.")
x <- 1 / sqrt(diag(Sigma))
R <- x * t(x * Sigma)}
else if(is.vector(Sigma)) {
k <- as.integer(sqrt(length(Sigma)))
Sigma <- matrix(Sigma, k, k)
x <- 1 / sqrt(diag(Sigma))
R <- as.vector(x * t(x * Sigma))}
return(R)
}
CovEstim <- function(Model, parm, Data, Method="Hessian")
{
if(Method == "Hessian") {
VarCov <- try(-as.inverse(Hessian(Model, parm, Data)),
silent=TRUE)
if(!inherits(VarCov, "try-error"))
diag(VarCov)[which(diag(VarCov) <= 0)] <- .Machine\$double.eps
else {
cat("\nWARNING: Failure to solve matrix inversion of ",
"Approx. Hessian.\n", sep="")
cat("NOTE: Identity matrix is supplied instead.\n")
VarCov <- diag(length(parm))}
}
else if(Method == "Identity") VarCov <- diag(length(parm))
else if(Method == "OPG") {
if(is.null(Data[["X"]])) stop("X is required in the data.")
y <- TRUE
if(is.null(Data[["y"]])) {
y <- FALSE
if(is.null(Data[["Y"]]))
stop("y or Y is required in the data.")}
if(y == TRUE) {
if(length(Data[["y"]]) != nrow(Data[["X"]]))
stop("length of y differs from rows in X.")
}
else {
if(nrow(Data[["Y"]]) != nrow(Data[["X"]]))
stop("The number of rows differs in y and X.")}
LIV <- length(parm)
VarCov <- matrix(0, LIV, LIV)
for (i in 1:nrow(Data[["X"]])) {
Data.temp <- Data
Data.temp\$X <- Data.temp\$X[i,,drop=FALSE]
if(y == TRUE) Data.temp\$y <- Data.temp\$y[i]
else Data.temp\$Y <- Data.temp\$Y[i,]
g <- partial(Model, parm, Data.temp)
VarCov <- VarCov + tcrossprod(g,g)}
VarCov <- as.inverse(as.symmetric.matrix(VarCov))
}
else if(Method == "Sandwich") {
B <- as.inverse(Hessian(Model, parm, Data))
if(is.null(Data[["X"]])) stop("X is required in the data.")
y <- TRUE
if(is.null(Data[["y"]])) {
y <- FALSE
if(is.null(Data[["Y"]]))
stop("y or Y is required in the data.")}
if(y == TRUE) {
if(length(Data[["y"]]) != nrow(Data[["X"]]))
stop("length of y differs from rows in X.")
}
else {
if(nrow(Data[["Y"]]) != nrow(Data[["X"]]))
stop("The number of rows differs in y and X.")}
LIV <- length(parm)
M <- matrix(0, LIV, LIV)
n <- nrow(Data[["X"]])
for (i in 1:n) {
Data.temp <- Data
Data.temp\$X <- Data.temp\$X[i,,drop=FALSE]
if(y == TRUE) Data.temp\$y <- Data.temp\$y[i]
else Data.temp\$Y <- Data.temp\$Y[i,]
g <- partial(Model, parm, Data.temp)
M <- M + tcrossprod(g,g)}
M <- as.symmetric.matrix(M)
}
else cat("\nWARNING: CovEst Method is unrecognized.")
return(VarCov)
}
GaussHermiteCubeRule <- function(N, dims, rule)
{
else Q <- rule
if(dims == 1) return(Q)
patterns_eq <- function(N, dims)
{
I <- matrix(1:N)
for (i in 2:dims) {
nf <- dim(I)[1]
nc <- dim(I)[2]
I2 <- cbind(kronecker(matrix(I[1, ], 1, nc),
matrix(1, N, 1)), (1:N))
for (j in 2:nf)
I2 <- rbind(I2, cbind(kronecker(matrix(I[j, ], 1, nc),
matrix(1, N, 1)), (1:N)))
I <- I2}
return(I)
}
I <- patterns_eq(N, dims)
n <- dim(I)[1]
X2 <- matrix(0, n, dims)
A2 <- matrix(1, n, 1)
for (i in 1:n)
for (j in 1:dims) {
X2[i, j] <- Q\$nodes[I[i, j]]
A2[i, 1] <- A2[i, 1] * Q\$weights[I[i, j]]}
Max <- Q\$weights[1] * Q\$weights[round((N + 1)/2)]/15
keep <- (A2 > Max)
n2 <- sum(keep)
X <- matrix(0, n2, dims)
A <- matrix(1, n2, 1)
k <- 0
for (i in 1:n)
if(keep[i]) {
k <- k + 1
X[k, ] <- X2[i, ]
A[k, ] <- A2[i, ]}
out <- list(nodes=X, weights=as.vector(A))
class(out) <- "gausshermitecuberule"
return(out)
}
Hessian <- function(Model, parm, Data, Interval=1e-6, Method="Richardson")
{
if(Method == "simple") {
parm.len <- length(parm)
eps <- Interval * parm
H <- matrix(0, parm.len, parm.len)
for (i in 1:parm.len) {
for (j in i:parm.len) {
x1 <- x2 <- x3 <- x4 <- parm
x1[i] <- x1[i] + eps[i]
x1[j] <- x1[j] + eps[j]
x2[i] <- x2[i] + eps[i]
x2[j] <- x2[j] - eps[j]
x3[i] <- x3[i] - eps[i]
x3[j] <- x3[j] + eps[j]
x4[i] <- x4[i] - eps[i]
x4[j] <- x4[j] - eps[j]
H[i, j] <- {Model(x1, Data)[["LP"]] -
Model(x2, Data)[["LP"]] - Model(x3, Data)[["LP"]] +
Model(x4, Data)[["LP"]]} / {4 * eps[i] * eps[j]}
}
}
H[lower.tri(H)] <- t(H)[lower.tri(H)]
return(H)
}
else if(Method != "Richardson") stop("Method is unknown.")
genD <- function(Model, parm, Data, Interval)
{
d <- 0.0001
r <- 4
v <- 2
zero.tol <- sqrt(.Machine\$double.eps / 7e-7)
f0 <- Model(parm, Data)[["LP"]]
p <- length(parm)
h0 <- abs(d*parm) + Interval*(abs(parm) < zero.tol)
D <- matrix(0, length(f0), (p*(p + 3)) / 2)
Daprox <- matrix(0, length(f0), r)
Hdiag  <- matrix(0, length(f0), p)
Haprox <- matrix(0, length(f0), r)
for (i in 1:p) {
h <- h0
for (k in 1:r) {
f1 <- Model(parm + (i == (1:p))*h, Data)[["LP"]]
f2 <- Model(parm - (i == (1:p))*h, Data)[["LP"]]
Daprox[,k] <- (f1 - f2)  / (2*h[i])
Haprox[,k] <- (f1 - 2*f0 + f2) / h[i]^2
h <- h / v
NULL}
for (m in 1:(r - 1))
for (k in 1:(r - m)) {
Daprox[,k] <- {Daprox[,k+1]*(4^m) -
Daprox[,k]} / (4^m - 1)
Haprox[,k] <- {Haprox[,k+1]*(4^m) -
Haprox[,k]} / (4^m - 1)
NULL}
D[,i] <- Daprox[,1]
Hdiag[,i] <- Haprox[,1]
NULL}
u <- p
for (i in 1:p) {
for (j in 1:i) {
u <- u + 1
if(i == j) {D[,u] <- Hdiag[,i]; NULL}
else {
h <- h0
for (k in 1:r) {
f1 <- Model(parm + (i == (1:p))*h +
(j == (1:p))*h, Data)[["LP"]]
f2 <- Model(parm - (i == (1:p))*h -
(j == (1:p))*h, Data)[["LP"]]
Daprox[,k] <- {f1 - 2*f0 + f2 -
Hdiag[,i]*h[i]^2 -
Hdiag[,j]*h[j]^2} / (2*h[i]*h[j])
h <- h / v}
for (m in 1:(r - 1))
for (k in 1:(r-m)) {
Daprox[,k] <- {Daprox[,k+1]*(4^m) -
Daprox[,k]} / (4^m - 1); NULL}
D[,u] <- Daprox[,1]
NULL}}}
invisible(D)
}
D <- genD(Model, parm, Data, Interval)
if(1 != nrow(D)) stop("BUG! should not get here.")
H <- diag(NA, length(parm))
u <- length(parm)
for (i in 1:length(parm)) {
for (j in 1:i) {
u <- u + 1
H[i,j] <- D[,u]}}
H <- H + t(H)
diag(H) <- diag(H) / 2
return(H)
}
is.positive.definite <- function(x)
{
if(!is.matrix(x)) stop("x is not a matrix.")
if(!is.square.matrix(x)) stop("x is not a square matrix.")
if(!is.symmetric.matrix(x)) stop("x is not a symmetric matrix.")
### Deprecated Method 1
#pd <- TRUE
#ed <- eigen(x, symmetric=TRUE)
#ev <- ed\$values
#if(!all(ev >= -1e-06 * abs(ev[1]))) pd <- FALSE
### Deprecated Method 2
#eval <- eigen(x, only.values=TRUE)\$values
#if(any(is.complex(eval))) eval <- rep(0, max(dim(x)))
#tol <- max(dim(x)) * max(abs(eval)) * .Machine\$double.eps
#if(all(eval > tol)) pd <- TRUE
#else pd <- FALSE
### Currently Active, Method 3
eigs <- eigen(x, symmetric=TRUE)\$values
if(any(is.complex(eigs))) return(FALSE)
if(all(eigs > 0)) pd <- TRUE
else pd <- FALSE
return(pd)
}
is.positive.semidefinite <- function(x)
{
if(!is.matrix(x)) stop("x is not a matrix.")
if(!is.square.matrix(x)) stop("x is not a square matrix.")
if(!is.symmetric.matrix(x)) stop("x is not a symmetric matrix.")
eigs <- eigen(x, symmetric=TRUE)\$values
if(any(is.complex(eigs))) return(FALSE)
if(all(eigs >= 0)) pd <- TRUE
else pd <- FALSE
return(pd)
}
is.square.matrix <- function(x) {return(nrow(x) == ncol(x))}
is.symmetric.matrix <- function(x) {return(sum(x == t(x)) == (nrow(x)^2))}
Jacobian <- function(Model, parm, Data, Interval=1e-6, Method="simple")
{
f <- Model(parm, Data)[[1]]
n <- length(parm)
if(Method == "simple") {
df <-matrix(NA, length(f), n)
for (i in 1:n) {
dx <- parm
dx[i] <- dx[i] + Interval
df[,i] <- {Model(dx, Data)[[1]] - f} / Interval}
return(df)
}
else if(Method == "Richardson") {
d <- 0.0001
zero.tol <- sqrt(.Machine\$double.eps / 7e-7)
r <- 4
v <- 2
a <- array(NA, c(length(f), r, n))
h <- abs(d*parm) + Interval*(abs(parm) < zero.tol)
for (k in 1:r) {
for (i in 1:n) {
a[,k,i] <- {Model(parm + h*(i == seq(n)), Data)[[1]] -
Model(parm - h*(i == seq(n)), Data)[[1]]} / (2*h[i])}
h <- h / v}
for (m in 1:(r - 1))
a <- {a[,2:(r+1-m),,drop=FALSE]*(4^m) -
a[,1:(r-m),,drop=FALSE]} / (4^m-1)
return(array(a, dim(a)[c(1,3)]))
}
else stop("The", Method, "is unknown.")
}
logdet <- function(x)
{
return(2*sum(log(diag(chol(x)))))
}
lower.triangle <- function(x, diag=FALSE)
{
return(x[lower.tri(x, diag=diag)])
}
size=0, na.rm=FALSE) {
if(nrow <= 0) nrow <- length(count.fields(file))
con <- file(file, open="r")
on.exit(close(con))
if(size <= 0) size <- nrow
if(samples <= 0) samples <- nrow
use <- sort(sample(nrow, samples))
ncol <- length(now)
col.names <- now
skip <- 1
}
else {
col.names <- paste("X[,", 1:ncol, "]", sep="")
skip <- 0
}
seek(con, 0)
X <- matrix(0, nrow=samples, ncol=ncol)
rownames(X) <- use
left <- nrow
got <- 1
while (left > 0) {
now <- matrix(scan(file=con, sep=sep, skip=skip, n=size*ncol,
quiet=TRUE), ncol=ncol, byrow=TRUE)
print(dim(now))
want <- (begin:end)[begin:end %in% use] - read
if(length(want) > 0) {
nowdat <- now[want,]
newgot <- got + length(want) - 1
X[got:newgot,] <- nowdat
got <- newgot + 1}
left <- left - size
}
colnames(X) <- col.names
if(na.rm == TRUE) {
num.mis <- sum(is.na(X))
if(num.mis > 0) {
cat("\n", num.mis, "missing value(s) found.")
cat("\n", sum(complete.cases(X)),
"row(s) found with missing values.")}}
return(X)
}
.rowVars <- function(X)
{
N <- ncol(X)
Y <- X - matrix(rowMeans(X), nrow(X), N)
Z <- rowMeans(Y*Y)*N/{N-1}
return(Z)
}
SparseGrid <- function(J, K)
{
### Initial Checks
J <- max(abs(round(J)), 1)
K <- max(abs(round(K)), 1)
type <- "GQN" #Gauss-Hermite
GQN <- function(level)
{
switch(level, {
n <- c(0)
w <- c(1)
}, {
n <- c(1)
w <- c(0.5)
}, {
n <- c(0, 1.73205080756888)
w <- c(0.666666666666667, 0.166666666666667)
}, {
n <- c(0.741963784302726, 2.3344142183389800)
w <- c(0.454124145231931, 0.0458758547680685)
}, {
n <- c(0, 1.35562617997427, 2.85697001387281)
w <- c(0.533333333333333, 0.222075922005613,
0.0112574113277207)
}, {
n <- c(0.616706590192594, 1.88917587775371,
3.32425743355212)
w <- c(0.408828469556029, 0.0886157460419145,
0.00255578440205624)
}, {
n <- c(0, 1.15440539473997, 2.36675941073454,
3.75043971772574)
w <- c(0.4571428571428580, 0.240123178605013000,
0.0307571239675865, 0.000548268855972219)
}, {
n <- c(0.539079811351375, 1.63651904243511,
2.802485861287540, 4.14454718612589)
w <- c(0.373012257679077, 0.117239907661759,
0.00963522012078826, 0.000112614538375368)
}, {
n <- c(0, 1.02325566378913, 2.07684797867783,
3.205429002856470, 4.51274586339978)
w <- c(0.406349206349207, 0.244097502894939,
0.049916406765218, 0.00278914132123177,
2.23458440077466e-05)
}, {
n <- c(0.484935707515498, 1.46598909439116,
2.484325841638950, 3.58182348355193,
4.859462828332310)
w <- c(0.3446423349320190, 0.13548370298026700,
0.0191115805007703, 0.00075807093431222,
4.31065263071831e-06)
}, {
n <- c(0, 0.928868997381064, 1.87603502015485,
2.86512316064364, 3.93616660712998,
5.18800122437487)
w <- c(0.369408369408370000, 0.24224029987397000,
0.066138746071057600, 0.00672028523553727,
0.000195671930271223, 8.1218497902149e-07)
}, {
n <- c(0.444403001944139, 1.34037519715162,
2.259464451000800, 3.22370982877010,
4.271825847932280, 5.50090170446775)
w <- c(0.3216643615128300, 0.14696704804533000,
0.0291166879123641, 0.00220338068753318,
4.83718492259061e-05, 1.49992716763716e-07)
}, {
n <- c(0, 0.85667949351945, 1.72541837958824,
2.62068997343221, 3.56344438028163,
4.59139844893652, 5.80016725238650)
w <- c(0.340992340992341, 0.237871522964136,
0.0791689558604501, 0.0117705605059965,
0.000681236350442926, 1.15265965273339e-05,
2.7226276428059e-08)
}, {
n <- c(0.412590457954602, 1.24268895548546,
2.088344745701940, 2.96303657983867,
3.886924575059770, 4.89693639734556,
6.08740954690129)
w <- c(0.3026346268130190, 0.15408333984251400,
0.0386501088242534, 0.00442891910694741,
0.000200339553760744, 2.66099134406763e-06,
4.86816125774839e-09)
}, {
n <- c(0, 0.799129068324548, 1.60671006902873,
2.43243682700976, 3.28908242439877,
4.19620771126902, 5.19009359130478,
6.36394788882984)
w <- c(0.3182595182595180, 0.2324622936097320,
0.0894177953998444, 0.0173657744921376,
0.00156735750354996, 5.64214640518902e-05,
5.9754195979206e-07, 8.58964989963318e-10)
}, {
n <- c(0.386760604500557, 1.16382910055496,
1.951980345716330, 2.76024504763070,
3.600873624171550, 4.49295530252001,
5.472225705949340, 6.63087819839313)
w <- c(0.286568521238012, 0.158338372750949,
0.0472847523540141, 0.00726693760118474,
0.00052598492657391, 1.53000321624873e-05,
1.30947321628682e-07, 1.49781472316183e-10)
}, {
n <- c(0, 0.751842600703896, 1.50988330779674,
2.28101944025299, 3.07379717532819,
3.90006571719801, 4.77853158962998,
5.74446007865941, 6.88912243989533)
w <- c(0.299538370126608, 0.226706308468979,
0.0974063711627181, 0.0230866570257112,
0.00285894606228465, 0.000168491431551339,
4.01267944797987e-06, 2.80801611793058e-08,
2.58431491937492e-11)
}, {
n <- c(0.365245755507698, 1.0983955180915,
1.83977992150865, 2.59583368891124,
3.37473653577809, 4.1880202316294,
5.05407268544274, 6.0077459113596,
7.13946484914648)
w <- c(0.2727832346542880, 0.160685303893513,
0.0548966324802227, 0.0105165177519414,
0.00106548479629165, 5.1798961441162e-05,
1.02155239763698e-06, 5.90548847883655e-09,
4.41658876935871e-12)
}, {
n <- c(0, 0.71208504404238, 1.42887667607837,
2.15550276131694, 2.89805127651575,
3.66441654745064, 4.46587262683103,
5.32053637733604, 6.26289115651325,
7.38257902403043)
w <- c(0.28377319275152100, 0.220941712199144000,
0.10360365727614400, 0.028666691030118500,
0.00450723542034204, 0.000378502109414268,
1.53511459546667e-05, 2.53222003209287e-07,
1.22037084844748e-09, 7.48283005405723e-13)
}, {
n <- c(0.346964157081356, 1.04294534880275,
1.74524732081413, 2.45866361117237,
3.18901481655339, 3.94396735065732,
4.73458133404606, 5.5787388058932,
6.51059015701366, 7.61904854167976)
w <- c(0.260793063449555, 0.161739333984,
0.061506372063976, 0.013997837447101,
0.00183010313108049, 0.000128826279961929,
4.40212109023086e-06, 6.12749025998296e-08,
2.48206236231518e-10, 1.25780067243793e-13)
}, {
n <- c(0, 0.678045692440644, 1.35976582321123,
2.04910246825716, 2.75059298105237,
3.46984669047538, 4.21434398168842,
4.99496394478203, 5.82938200730447,
6.75144471871746, 7.84938289511382)
w <- c(0.270260183572877, 0.21533371569506,
0.108392285626419, 0.0339527297865428,
0.00643969705140878, 0.000708047795481537,
4.21923474255159e-05, 1.22535483614825e-06,
1.45066128449307e-08, 4.97536860412175e-11,
2.09899121956567e-14)
}, {
n <- c(0.331179315715274, 0.995162422271216,
1.66412483911791, 2.34175999628771,
3.03240422783168, 3.74149635026652,
4.47636197731087, 5.24772443371443,
6.07307495112290, 6.98598042401882,
8.07402998402171)
w <- c(0.250243596586935, 0.161906293413675,
0.0671963114288899, 0.0175690728808058,
0.00280876104757721, 0.000262283303255964,
1.33459771268087e-05, 3.319853749814e-07,
3.36651415945821e-09, 9.84137898234601e-12,
3.47946064787714e-15)
}, {
n <- c(0, 0.648471153534496, 1.29987646830398,
1.95732755293342, 2.62432363405918,
3.30504002175297, 4.00477532173330,
4.73072419745147, 5.49347398647179,
6.31034985444840, 7.21465943505186,
8.29338602741735)
w <- c(0.258509740808839, 0.209959669577543,
0.112073382602621, 0.0388671837034809,
0.00857967839146566, 0.00116762863749786,
9.3408186090313e-05, 4.08997724499215e-06,
8.77506248386172e-08, 7.67088886239991e-10,
1.92293531156779e-12, 5.73238316780209e-16)
}, {
n <- c(0.317370096629452, 0.953421922932109,
1.593480429816420, 2.240467851691750,
2.89772864322331, 3.56930676407356,
4.26038360501991, 4.97804137463912,
5.73274717525120, 6.54167500509863,
7.43789066602166, 8.50780351919526)
w <- c(0.240870115546641, 0.161459512867,
0.0720693640171784, 0.021126344408967,
0.00397660892918131, 0.000464718718779398,
3.2095005652746e-05, 1.21765974544258e-06,
2.26746167348047e-08, 1.71866492796487e-10,
3.71497415276242e-13, 9.39019368904192e-17)
}, {
n <- c(0, 0.622462279186076, 1.24731197561679,
1.87705836994784, 2.51447330395221,
3.16277567938819, 3.82590056997249,
4.50892992296729, 5.21884809364428,
5.96601469060670, 6.76746496380972,
7.65603795539308, 8.71759767839959)
w <- c(0.248169351176485, 0.20485102565034,
0.114880924303952, 0.043379970167645,
0.0108567559914623, 0.0017578504052638,
0.000177766906926527, 1.06721949052025e-05,
3.5301525602455e-07, 5.73802386889938e-09,
3.79115000047719e-11, 7.10210303700393e-14,
1.53003899799868e-17)
})
return(list(nodes=n, weights=w))
}
SparseGridGetSeq <- function(J, norm)
{
seq.vec <- rep(0, J)
a <- norm - J
seq.vec[1] <- a
fs <- matrix(seq.vec, nrow=1, ncol=length(seq.vec))
cnt <- 1
while (seq.vec[J] < a) {
if(cnt == J) {
for (i in seq(cnt - 1, 1, -1)) {
cnt <- i
if(seq.vec[i] != 0) {
break
}
}
}
seq.vec[cnt] <- seq.vec[cnt] - 1
cnt <- cnt + 1
seq.vec[cnt] <- a - sum(seq.vec[1:(cnt - 1)])
if(cnt < J) {
seq.vec[(cnt + 1):J] <- rep(0, J -
cnt)
}
fs <- rbind(fs, seq.vec)
}
fs <- fs + 1
return(fs)
}
SparseGridKronProd <- function(n1D, w1D)
{
nodes <- matrix(n1D[[1]], nrow=length(n1D[[1]]), ncol=1)
weights <- w1D[[1]]
if(length(n1D) > 1) {
for (j in 2:length(n1D)) {
newnodes <- n1D[[j]]
nodes <- cbind(kronecker(nodes, rep(1, length(newnodes))),
kronecker(rep(1, nrow(nodes)), newnodes))
weights <- kronecker(weights, w1D[[j]])
}
}
return(list(nodes=nodes, weights=weights))
}
sortrows <- function (A, index.return=FALSE)
{
if(!is.matrix(A)) {
stop("SparseGrid:::sortrows expects a matrix as argument A.")
}
A.nrow <- nrow(A)
A.ncol <- ncol(A)
if(index.return == TRUE) indices <- 1:nrow(A)
for (col.cnt in seq(ncol(A), 1, -1)) {
tmp.indices <- order(A[, col.cnt])
A <- A[tmp.indices, , drop=FALSE]
if(index.return == TRUE) indices <- indices[tmp.indices]
}
if(index.return == TRUE) {
res <- list(x=matrix(A, nrow=A.nrow, ncol=A.ncol),
ix=indices)
}
else res <- matrix(A, nrow=A.nrow, ncol=A.ncol)
return(res)
}
tryCatch({
n1D <- vector(mode="list", length=K)
w1D <- vector(mode="list", length=K)
R1D <- rep(0, K)
for (level in 1:K) {
res <- eval(call(type, level))
nodes <- res\$nodes
weights <- res\$weights
R1D[level] <- length(weights)
n1D[[level]] <- nodes
w1D[[level]] <- weights
}
}, error=function(e) {cat("Error evaluating the 1D rule\n")})
minq <- max(0, K - J)
maxq <- K - 1
nodes <- matrix(0, nrow=0, ncol=J)
weights <- numeric(0)
for (q.cnt in minq:maxq) {
r <- length(weights)
bq <- (-1)^(maxq - q.cnt) * choose(J - 1, J +
q.cnt - K)
indices.mat <- SparseGridGetSeq(J, J + q.cnt)
Rq <- sapply(1:nrow(indices.mat), function(row.cnt) {
prod(R1D[indices.mat[row.cnt, ]])})
Rq.sum <- sum(Rq)
nodes <- rbind(nodes, matrix(0, nrow=Rq.sum, ncol=J))
weights <- c(weights, rep(0, Rq.sum))
for (j in 1:nrow(indices.mat)) {
midx <- indices.mat[j, ]
res <- SparseGridKronProd(n1D[midx], w1D[midx])
nodes[(r + 1):(r + Rq[j]), ] <- res\$nodes
weights[(r + 1):(r + Rq[j])] <- bq * res\$weights
r <- r + Rq[j]
}
nodes.sorted <- sortrows(nodes, index.return=TRUE)
nodes <- nodes.sorted\$x
weights <- weights[nodes.sorted\$ix]
keep <- 1
lastkeep <- 1
if(nrow(nodes) > 1) {
for (j in 2:nrow(nodes)) {
if(all(nodes[j, ] == nodes[j - 1, ])) {
weights[lastkeep] <- weights[lastkeep] + weights[j]
}
else {
lastkeep <- j
keep <- c(keep, j)
}
}
}
nodes <- matrix(nodes[keep, ], nrow=length(keep), ncol=J)
weights <- weights[keep]
}
nr <- length(weights)
m <- n1D[[1]]
for (j in 1:J) {
keep <- rep(0, nr)
numnew <- 0
for (r in 1:nr) {
if(nodes[r, j] != m) {
numnew <- numnew + 1
keep[numnew] <- r
}
}
if(numnew > 0) {
nodes <- rbind(nodes, matrix(nodes[keep[1:numnew],
], nrow=numnew, ncol=J))
nodes[(nr + 1):(nr + numnew), j] <- 2 * m - nodes[(nr +
1):(nr + numnew), j]
weights <- c(weights, weights[keep[1:numnew]])
nr <- nr + numnew
}
}
nodes.sorted <- sortrows(nodes, index.return=TRUE)
nodes <- nodes.sorted\$x
weights <- weights[nodes.sorted\$ix]
weights <- abs(weights)/sum(abs(weights))
return(list(nodes=nodes, weights=weights))
}
TransitionMatrix <- function(theta.y=NULL, y.theta=NULL, p.theta=NULL)
{
if(!is.null(theta.y)) {
theta.y <- as.vector(theta.y)
N <- length(theta.y)
theta.y <- as.matrix(table(theta.y[-N], theta.y[-1]))
theta.y <- theta.y / rowSums(theta.y)
return(theta.y)
}
else {
N <- length(y.theta)
y.theta <- as.matrix(table(From=y.theta[-N], To=y.theta[-1]))
y.theta <- y.theta / rowSums(y.theta)
if(!is.null(p.theta)) {
if(!is.matrix(p.theta)) p.theta <- as.matrix(p.theta)
if(!identical(dim(y.theta),dim(p.theta)))
stop("Dimensions of p.theta differ from the matrix of y.theta.")
theta.y <- y.theta * p.theta
theta.y <- theta.y / rowSums(theta.y)
}
else p.theta <- y.theta
return(p.theta)
}
}
tr <- function(x) {return(sum(diag(x)))}
upper.triangle <- function(x, diag=FALSE)
{
return(x[upper.tri(x, diag=diag)])
}

#End
```

## Try the LaplacesDemon package in your browser

Any scripts or data that you put into this service are public.

LaplacesDemon documentation built on July 1, 2018, 9:02 a.m.