Description Usage Arguments Value Author(s) References Examples
See above
1 2 | generalizedJointF(x, Vx, mu, s, rho_X, G_Point7,
GH_Quadrature, maxSize, choleskySpeed = TRUE, cores = 1, verbose = -1)
|
x |
List of vectors of points to estimate the fitted numerical density functions at. |
Vx |
List of samples from each vector |
mu |
List of means of each vector |
s |
List of standard deviations of each vector |
rho_X |
Correlation matrix of the vectors |
G_Point7 |
Seven Points of Gaussian Hermite Quadrature |
GH_Quadrature |
Seven Weights of Gaussian Hermite Quadrature |
maxSize |
Used for a speed improvement, the larger this is the faster the function will work as it iterates over each bin, but within each bin vectorized operations take place. The restriction is memory. Typically floor(sqrt(1450000)) performs well on 32 Bit machines, but it is up to user discretion to find the optimal value for their machine. |
choleskySpeed |
TRUE to apply a slight speed improvement via a cholesky decomposition. Occasionally this will not work if the matrices are not able to be decomposed in such a manner, theoretically impossible for a correlation matrix but can occur with matrices that are almost singular. |
cores |
Number of cores to use in parallel processing, only works with 1 core on windows machines. |
verbose |
If greater than 1 it provides additional information to the console. |
An array of dimension corresponding to the dimension of each vector in x.
Luke Mazur
no references
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 | ##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function (x, Vx, mu, s, rho_X, G_Point7, GH_Quadrature, maxSize,
choleskySpeed = TRUE, cores = 1, verbose = -1)
{
numberOfVectors <- length(x)
if (length(Vx) != numberOfVectors) {
stop("x and Vx must be the same length")
}
numberOfVectors <- length(x)
if (length(mu) != numberOfVectors) {
stop("x and mu must be the same length")
}
numberOfVectors <- length(x)
if (length(s) != numberOfVectors) {
stop("x and s must be the same length")
}
n <- unlist(lapply(1:numberOfVectors, FUN = function(i) {
return(length(x[[i]]))
}))
fhat <- lapply(1:numberOfVectors, FUN = function(i) {
return(kde(x = Vx[[i]], binned = TRUE))
})
if (verbose > 1) {
print("calculating Nataf_rho matrix")
}
Nataf_rho <- matrix(rep(1, (numberOfVectors^2)), nrow = numberOfVectors,
ncol = numberOfVectors)
Nataf_rho[upper.tri(Nataf_rho, diag = TRUE)] <- NA
Nataf_rho <- mclapply(1:numberOfVectors, mc.cores = cores,
FUN = function(j) {
return(lapply(1:numberOfVectors, FUN = function(i) {
if (verbose > 1) {
print("row and column")
print(i)
print(j)
}
if (!is.na(Nataf_rho[i, j])) {
return(rho_0(Vx[[j]], Vx[[i]], mu[[j]], mu[[i]],
s[[j]], s[[i]], rho_X[i, j], G_Point7, GH_Quadrature,
verbose))
} else {
return(NA)
}
}))
})
Nataf_rho <- unlist(Nataf_rho)
Nataf_rho <- matrix(Nataf_rho, nrow = numberOfVectors, ncol = numberOfVectors)
diag(Nataf_rho) <- 1
Nataf_rho[upper.tri(Nataf_rho, diag = FALSE)] <- Nataf_rho[lower.tri(Nataf_rho,
diag = FALSE)]
if (verbose > 1) {
print("calculating y and phiy")
}
y <- lapply(1:numberOfVectors, FUN = function(i) {
return(qnorm(pkde(x[[i]], fhat[[i]])))
})
phiy <- lapply(1:numberOfVectors, FUN = function(i) {
return(exp(-((y[[i]])^2)/2)/sqrt(2 * pi))
})
if (verbose > 1) {
print("calculating inverse of Nataf_rho")
}
A <- solve(Nataf_rho)
if (verbose > 1) {
print("calculating all possible combinations")
}
allPossibleCombinationsX <- expand.grid(x)
allPossibleCombinationsY <- expand.grid(y)
allPossibleCombinationsPhiy <- expand.grid(phiy)
yIndependentPartOfPhi <- (1/(sqrt((2 * pi)^numberOfVectors *
det(Nataf_rho))))
nrowAllPossibleCombinationsY <- nrow(allPossibleCombinationsY)
numberOfBins <- ceiling(nrowAllPossibleCombinationsY/maxSize)
binFactors <- cut(1:nrowAllPossibleCombinationsY, breaks = numberOfBins)
bins <- split(allPossibleCombinationsY, f = binFactors)
if (verbose > 1) {
print(paste("There are ", numberOfBins, " bins."))
}
if (choleskySpeed) {
cholA <- chol(A)
yDependentPartOfPhi <- unlist(mclapply(1:numberOfBins,
mc.cores = cores, FUN = function(i) {
if (verbose > 1) {
print(i)
}
temp <- as.matrix(bins[[i]])
temp <- crossprod(cholA %*% t(temp))
return(exp(-diag(temp)))
}))
}
else {
yDependentPartOfPhi <- unlist(mclapply(1:numberOfBins,
mc.cores = cores, FUN = function(i) {
if (verbose > 1) {
print(i)
}
temp <- as.matrix(bins[[i]])
temp <- temp %*% A %*% t(temp)
return(exp(-diag(temp)))
}))
}
if (verbose > 1) {
print("Moving on to xAndPhiy")
}
xAndPhiy <- mclapply(1:numberOfVectors, mc.cores = cores,
FUN = function(j) {
return(dkde(allPossibleCombinationsX[, j], fhat[[j]])/allPossibleCombinationsPhiy[,
j])
})
if (verbose > 1) {
print("Taking products of xAndPhiy")
}
xAndPhiy <- unlist(apply(as.data.frame(xAndPhiy), 1, FUN = prod))
f <- yIndependentPartOfPhi * yDependentPartOfPhi * xAndPhiy
f <- array(f, dim = n)
return(f)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.