# Consider multivariate normal vector:
# X = (X1, X2, X3, X4, X5) ~ N(mean, sigma)
# Prepare multivariate normal vector parameters
# expected value
mean <- c(-2, -1, 0, 1, 2)
n_dim <- length(mean)
# correlation matrix
cor <- c( 1, 0.1, 0.2, 0.3, 0.4,
0.1, 1, -0.1, -0.2, -0.3,
0.2, -0.1, 1, 0.3, 0.2,
0.3, -0.2, 0.3, 1, -0.05,
0.4, -0.3, 0.2, -0.05, 1)
cor <- matrix(cor, ncol = n_dim, nrow = n_dim, byrow = TRUE)
# covariance matrix
sd_mat <- diag(c(1, 1.5, 2, 2.5, 3))
sigma <- sd_mat %*% cor %*% t(sd_mat)
# Estimate probability:
# P(-3 < X1 < 1, -2.5 < X2 < 1.5, -2 < X3 < 2, -1.5 < X4 < 2.5, -1 < X5 < 3)
lower <- c(-3, -2.5, -2, -1.5, -1)
upper <- c(1, 1.5, 2, 2.5, 3)
p.list <- pmnorm(lower = lower, upper = upper,
mean = mean, sigma = sigma)
p <- p.list$prob
print(p)
# Additionally estimate a probability
lower.1 <- c(-Inf, 0, -Inf, 1, -Inf)
upper.1 <- c(Inf, Inf, 3, 4, 5)
lower.mat <- rbind(lower, lower.1)
upper.mat <- rbind(upper, upper.1)
p.list.1 <- pmnorm(lower = lower.mat, upper = upper.mat,
mean = mean, sigma = sigma)
p.1 <- p.list.1$prob
print(p.1)
# Estimate the probabilities
# P(-1 < X1 < 1, -3 < X3 < 3, -5 < X5 < 5 | X2 = -2, X4 = 4)
lower.2 <- c(-1, -3, -5)
upper.2 <- c(1, 3, 5)
given_ind <- c(2, 4)
given_x <- c(-2, 4)
p.list.2 <- pmnorm(lower = lower.2, upper = upper.2,
mean = mean, sigma = sigma,
given_ind = given_ind, given_x = given_x)
p.2 <- p.list.2$prob
print(p.2)
# Additionally estimate the probability
# P(-Inf < X1 < 1, -3 < X3 < Inf, -Inf < X5 < Inf | X2 = 4, X4 = -2)
lower.3 <- c(-Inf, -3, -Inf)
upper.3 <- c(1, Inf, Inf)
given_x.1 <- c(-2, 4)
lower.mat.2 <- rbind(lower.2, lower.3)
upper.mat.2 <- rbind(upper.2, upper.3)
given_x.mat <- rbind(given_x, given_x.1)
p.list.3 <- pmnorm(lower = lower.mat.2, upper = upper.mat.2,
mean = mean, sigma = sigma,
given_ind = given_ind, given_x = given_x.mat)
p.3 <- p.list.3$prob
print(p.3)
# Estimate the gradient of previous probabilities respect various arguments
p.list.4 <- pmnorm(lower = lower.mat.2, upper = upper.mat.2,
mean = mean, sigma = sigma,
given_ind = given_ind, given_x = given_x.mat,
grad_lower = TRUE, grad_upper = TRUE,
grad_sigma = TRUE, grad_given = TRUE)
p.4 <- p.list.4$prob
print(p.4)
# Gradient respect to 'lower'
grad_lower <- p.list.4$grad_lower
# for the first probability
print(grad_lower[1, ])
# for the second probability
print(grad_lower[2, ])
# Gradient respect to 'upper'
grad_upper <- p.list.4$grad_upper
# for the first probability
print(grad_upper[1, ])
# for the second probability
print(grad_upper[2, ])
# Gradient respect to 'given_x'
grad_given <- p.list.4$grad_given
# for the first probability
print(grad_given[1, ])
# for the second probability
print(grad_given[2, ])
# Gradient respect to 'sigma'
grad_given <- p.list.4$grad_given
# for the first probability
print(grad_given[1, ])
# for the second probability
print(grad_given[2, ])
# Compare analytical gradients from the previous example with
# their numeric (forward difference) analogues for the first probability
n_dependent <- length(lower.2)
n_given <- length(given_x)
n_dim <- n_dependent + n_given
delta <- 1e-6
grad_lower.num <- rep(NA, n_dependent)
grad_upper.num <- rep(NA, n_dependent)
grad_given.num <- rep(NA, n_given)
grad_sigma.num <- matrix(NA, nrow = n_dim, ncol = n_dim)
for (i in 1:n_dependent)
{
# respect to lower
lower.delta <- lower.2
lower.delta[i] <- lower.2[i] + delta
p.list.delta <- pmnorm(lower = lower.delta, upper = upper.2,
given_x = given_x,
mean = mean, sigma = sigma,
given_ind = given_ind)
grad_lower.num[i] <- (p.list.delta$prob - p.list.4$prob[1]) / delta
# respect to upper
upper.delta <- upper.2
upper.delta[i] <- upper.2[i] + delta
p.list.delta <- pmnorm(lower = lower.2, upper = upper.delta,
given_x = given_x,
mean = mean, sigma = sigma,
given_ind = given_ind)
grad_upper.num[i] <- (p.list.delta$prob - p.list.4$prob[1]) / delta
}
for (i in 1:n_given)
{
# respect to lower
given_x.delta <- given_x
given_x.delta[i] <- given_x[i] + delta
p.list.delta <- pmnorm(lower = lower.2, upper = upper.2,
given_x = given_x.delta,
mean = mean, sigma = sigma,
given_ind = given_ind)
grad_given.num[i] <- (p.list.delta$prob - p.list.4$prob[1]) / delta
}
for (i in 1:n_dim)
{
for(j in 1:n_dim)
{
# respect to sigma
sigma.delta <- sigma
sigma.delta[i, j] <- sigma[i, j] + delta
sigma.delta[j, i] <- sigma[j, i] + delta
p.list.delta <- pmnorm(lower = lower.2, upper = upper.2,
given_x = given_x,
mean = mean, sigma = sigma.delta,
given_ind = given_ind)
grad_sigma.num[i, j] <- (p.list.delta$prob - p.list.4$prob[1]) / delta
}
}
# Comparison of gradients respect to lower integration limits
h.lower <- cbind(analytical = p.list.4$grad_lower[1, ],
numeric = grad_lower.num)
print(h.lower)
# Comparison of gradients respect to upper integration limits
h.upper <- cbind(analytical = p.list.4$grad_upper[1, ],
numeric = grad_upper.num)
print(h.upper)
# Comparison of gradients respect to given values
h.given <- cbind(analytical = p.list.4$grad_given[1, ],
numeric = grad_given.num)
print(h.given)
# Comparison of gradients respect to the covariance matrix
h.sigma <- list(analytical = p.list.4$grad_sigma[, , 1],
numeric = grad_sigma.num)
print(h.sigma)
# Let's again estimate probability
# P(-1 < X1 < 1, -3 < X3 < 3, -5 < X5 < 5 | X2 = -2, X4 = 4)
# But assume that standardized (to zero mean and unit variance):
# 1) X1 and X2 have standardized PGN distribution with coefficients vectors
# pc1 = (0.5, -0.2, 0.1) and pc2 = (0.2, 0.05) correspondingly.
# 2) X3 has standardized student distribution with 5 degrees of freedom
# 3) X4 has standardized logistic distribution
# 4) X5 has standard normal distribution
marginal <- list(PGN = c(0.5, -0.2, 0.1), hpa = c(0.2, 0.05),
student = 5, logistic = numeric(), normal = NULL)
p.list.5 <- pmnorm(lower = lower.2, upper = upper.2,
mean = mean, sigma = sigma,
given_ind = given_ind, given_x = given_x,
grad_lower = TRUE, grad_upper = TRUE,
grad_sigma = TRUE, grad_given = TRUE,
marginal = marginal, grad_marginal = TRUE)
# Lets investigate derivatives respect to parameters
# of marginal distributions
# respect to pc1 of X1
p.list.5$grad_marginal[[1]]
# respect to pc2 of X2
p.list.5$grad_marginal[[2]]
# derivative respect to degrees of freedom of X5 is
# currently unavailable and will be set to 0
p.list.5$grad_marginal[[3]]
Run the code above in your browser using DataLab