# R/ic.weights.R In ic.infer: Inequality Constrained Inference in Linear Normal Situations

```ic.weights <- function (corr, ...)
{
## ... given to pmvnorm, should contain algorithm information only
## changed to ... in order to stay tuned with updates in mvtnorm
## and be also compatible with old versions (e.g. 2.5.0)
if (!(is.matrix(corr) & nrow(corr) == ncol(corr)))
stop("corr must be a square matrix.")
if (!(all(eigen(corr)\$value > 0)))
stop("corr must be positive definite.")
g <- nrow(corr)
liste <- 1:g
## weights in the order of highest to lowest dimension
## for variable with correlation or cov corr
weights <- rep(0, g + 1)
names(weights) <- rev(0:g)
weights[1] <- pmvnorm(lower = rep(0, g), upper = rep(Inf,
g), sigma = corr, ...)
weights[g + 1] <- pmvnorm(lower = rep(0, g), upper = rep(Inf,
g), sigma = solve(corr), ...)
## prevent lengthy calculations
## if longest vector too long for available storage
if (g > 4){
if (!is.numeric(try(matrix(0, floor((g-2)/2), choose(g, floor((g-2)/2))),
silent = TRUE)))
stop(paste("ic.weights will not work, corr too large, \n",
"interim matrix with ", floor((g-2)/2) * choose(g, floor((g-2)/2)),
" elements cannot be created.", sep = ""))
}
if (g==2) weights[2] <- 1-sum(weights)
if (g==3) {
weights[2] <- 0.5 - weights[4]
weights[3] <- 0.5 - weights[1]
}
if (g > 3) {
for (k in 1:floor((g - 2)/2)) {
### fill weights simultaneously from top and bottom
### fill the two middle weights by 0.5 - odd and even sum of others,
###     respectively, according to Silvapulle and Sen 2004,
###     Prop. 3.6.1, part 3
jetzt <- nchoosek(g, k)
wjetzt <- matrix(0, choose(g, k), 2)
for (j in 1:(choose(g, k))) {
diese <- jetzt[, j]
andere <- setdiff(liste, diese)
hilf <- corr[andere, andere] - corr[andere, diese] %*%
solve(corr[diese, diese], matrix(corr[diese,
andere], k, g - k))
## matrix necessary, because vector otherwise wrong direction
wjetzt[j, 1] <- pmvnorm(lower = rep(0, k), upper = rep(Inf,
k), sigma = solve(corr[diese, diese]), ...) *
pmvnorm(lower = rep(0, g - k), upper = rep(Inf,
g - k), sigma = hilf, ...)
if (!k==(g-2)/2){
## needed for odd g only for the last one from top,
## for even g calculated as 0.5 minus the other even weights
hilf <- corr[diese, diese] - corr[diese, andere] %*%
solve(corr[andere, andere], matrix(corr[andere,
diese], g - k, k))
wjetzt[j, 2] <- pmvnorm(lower = rep(0, g - k),
upper = rep(Inf, g - k), sigma = solve(corr[andere,
andere]), ...) * pmvnorm(lower = rep(0, k),
upper = rep(Inf, k), sigma = hilf, ...)
}
}
weights[k + 1] <- sum(wjetzt[, 1])
weights[g + 1 - k] <- sum(wjetzt[, 2])
}
### fill last weight using sum of odd and even weights
### even weights in odd positions and vice versa
if (g/2 == floor(g/2)) {
even.sum <- sum(weights[1+2*((g/2):0)])
odd.sum <- sum(weights[2*((g/2):1)])
if (g/4 == floor(g/4)) {
weights[g/2 + 1] <- 0.5 - even.sum ## even weights
weights[g/2 + 2] <- 0.5 - odd.sum ## odd weights
}
else {
weights[g/2 + 1] <- 0.5 - odd.sum ## odd weights
weights[g/2 + 2] <- 0.5 - even.sum ## even weights
}
}
else {
even.sum <- sum(weights[2*(((g+1)/2):1)])
odd.sum <- sum(weights[2*(((g+1)/2):1)-1])
if ((g+1)/4 == floor((g+1)/4)) {
weights[(g + 1)/2] <- 0.5 - even.sum ## even weights
weights[(g + 3)/2] <- 0.5 - odd.sum ## odd weights
}
else {
weights[(g + 1)/2] <- 0.5 - odd.sum ## odd weights
weights[(g + 3)/2] <- 0.5 - even.sum ## even weights
}
}
}
weights
}
```

## Try the ic.infer package in your browser

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

ic.infer documentation built on Oct. 5, 2023, 5:09 p.m.