[Vinecopula-commits] r82 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fr Feb 20 14:51:02 CET 2015
Author: tnagler
Date: 2015-02-20 14:51:01 +0100 (Fri, 20 Feb 2015)
New Revision: 82
Modified:
pkg/R/BiCopEst.r
pkg/R/RVineCopSelect.r
pkg/R/RVineSim.R
Log:
Bugfixes:
- BiCopEst: extend search interval for Tawn MLE to avoid optim-errors
- RVineSim: reorder U so that it corresponds to the order of RVM
Modified: pkg/R/BiCopEst.r
===================================================================
--- pkg/R/BiCopEst.r 2015-02-20 11:46:26 UTC (rev 81)
+++ pkg/R/BiCopEst.r 2015-02-20 13:51:01 UTC (rev 82)
@@ -898,17 +898,18 @@
MLE_intern_Tawn <- function(data, start.parm, family, se = FALSE) {
n <- dim(data)[1]
+
+ ## set bounds for optimization
tau <- fasttau(data[, 1], data[, 2])
-
if (family == 104 || family == 114 || family == 204 || family == 214) {
- parlower <- c(1.001, max(tau, 1e-04))
- parupper <- c(20, min(tau + 0.1, 0.99))
+ parlower <- c(1.001, max(tau - 0.1, 1e-04))
+ parupper <- c(20, min(tau + 0.2, 0.99))
} else if (family == 124 || family == 134 || family == 224 || family == 234) {
- parlower <- c(-20, max(-tau, 1e-04))
- parupper <- c(-1.001, min(-tau + 0.1, 0.99))
+ parlower <- c(-20, max(-tau - 0.1, 1e-04))
+ parupper <- c(-1.001, min(-tau + 0.2, 0.99))
}
- # Hier fehlt noch die log-likelihood Funktion
+ ## log-liklihood function
loglikfunc <- function(param) {
ll <- .C("LL_mod2",
as.integer(family),
@@ -925,6 +926,7 @@
return(ll)
}
+ ## optimize log-likelihood
out <- list()
# print(start.parm)
if (se == TRUE) {
@@ -951,6 +953,7 @@
control = list(fnscale = -1, maxit = 500))
}
+ ## return results
out$par <- optimout$par
out$value <- optimout$value
return(out)
Modified: pkg/R/RVineCopSelect.r
===================================================================
--- pkg/R/RVineCopSelect.r 2015-02-20 11:46:26 UTC (rev 81)
+++ pkg/R/RVineCopSelect.r 2015-02-20 13:51:01 UTC (rev 82)
@@ -1,13 +1,12 @@
RVineCopSelect <- function(data, familyset = NA, Matrix, selectioncrit = "AIC", indeptest = FALSE, level = 0.05, trunclevel = NA) {
-
n <- dim(data)[2]
N <- nrow(data)
+ ## sanity checks
if (dim(Matrix)[1] != dim(Matrix)[2])
stop("Structure matrix has to be quadratic.")
if (max(Matrix) > dim(Matrix)[1])
stop("Error in the structure matrix.")
-
if (N < 2)
stop("Number of observations has to be at least 2.")
if (n < 2)
@@ -24,59 +23,60 @@
if (level < 0 & level > 1)
stop("Significance level has to be between 0 and 1.")
+ ## adjustement for truncated vines
if (is.na(trunclevel))
trunclevel <- n
-
types <- familyset
if (trunclevel == 0)
types <- 0
+ ## reorder matrix to natural order
M <- Matrix
-
Mold <- M
-
o <- diag(M)
M <- reorderRVineMatrix(M)
-
data <- data[, o[length(o):1]]
+ ## create matrices required for selection of h-functions
MaxMat <- createMaxMat(M)
CondDistr <- neededCondDistr(M)
+ ## create objects for results
Types <- matrix(0, n, n)
-
Params <- matrix(0, n, n)
Params2 <- matrix(0, n, n)
-
V <- list()
V$direct <- array(NA, dim = c(n, n, N))
V$indirect <- array(NA, dim = c(n, n, N))
V$direct[n, , ] <- t(data[, n:1])
+ ## loop over all trees and pair-copulas
for (i in (n - 1):1) {
for (k in n:(i + 1)) {
- m <- MaxMat[k, i]
-
+ ## get (pseudo-) observations
+ m <- MaxMat[k, i] # edge identifier
zr1 <- V$direct[k, i, ]
-
if (m == M[k, i]) {
zr2 <- V$direct[k, (n - m + 1), ]
} else {
zr2 <- V$indirect[k, (n - m + 1), ]
}
+ ## estimate pair-copula
if (n + 1 - k > trunclevel) {
outcop <- BiCopSelect(zr2, zr1, 0, selectioncrit, indeptest, level)
} else {
# outcop = BiCopSelect(zr1,zr2,types,selectioncrit,indeptest,level)
- outcop <- BiCopSelect(zr2, zr1, types, selectioncrit, indeptest, level)
+ outcop <- BiCopSelect(zr2, zr1, types[1:3], selectioncrit, indeptest, level)
}
+ ## store results for pair-copula
Types[k, i] <- outcop$family
Params[k, i] <- outcop$par
Params2[k, i] <- outcop$par2
+ ## calculate pseudo observations required in the next tree
if (CondDistr$direct[k - 1, i])
# V$direct[k-1,i,] = outcop$CondOn.2
V$direct[k - 1, i, ] <- .C("Hfunc1",
@@ -102,15 +102,13 @@
}
}
+ ## return results
varnames <- paste("V", 1:n, sep = "")
-
print(Types)
-
RVM <- RVineMatrix(Mold,
family = Types,
par = Params,
par2 = Params2,
names = varnames)
-
return(RVM)
}
Modified: pkg/R/RVineSim.R
===================================================================
--- pkg/R/RVineSim.R 2015-02-20 11:46:26 UTC (rev 81)
+++ pkg/R/RVineSim.R 2015-02-20 13:51:01 UTC (rev 82)
@@ -1,21 +1,24 @@
RVineSim <- function(N, RVM, U = NULL) {
+
+ ## sanity checks
stopifnot(N >= 1)
if (!is(RVM, "RVineMatrix"))
stop("'RVM' has to be an RVineMatrix object.")
+ ## reorder matrix and U (if provided)
n <- dim(RVM)
-
o <- diag(RVM$Matrix)
RVM <- normalizeRVineMatrix(RVM)
-
takeU <- !is.null(U)
if (takeU) {
if (!is.matrix(U))
U <- rbind(U, deparse.level = 0L)
if ((d <- ncol(U)) < 2)
stop("U should be at least bivariate") # should be an (N, n) matrix
+ U <- U[, rev(o)]
}
+ ## create objects for C-call
matri <- as.vector(RVM$Matrix)
w1 <- as.vector(RVM$family)
th <- as.vector(RVM$par)
@@ -28,9 +31,9 @@
th2[is.na(th2)] <- 0
maxmat[is.na(maxmat)] <- 0
conindirect[is.na(conindirect)] <- 0
-
tmp <- rep(0, n * N)
+ ## simulate R-Vine
tmp <- .C("SimulateRVine",
as.integer(N),
as.integer(n),
@@ -45,13 +48,11 @@
as.integer(takeU),
PACKAGE = "VineCopula")[[9]]
+ ## store results, bring back to initial order and return
out <- matrix(tmp, ncol = n)
-
-
if (!is.null(RVM$names)) {
colnames(out) <- RVM$names
}
-
out <- out[, sort(o[length(o):1], index.return = TRUE)$ix]
return(out)
}
Mehr Informationen über die Mailingliste Vinecopula-commits