Capítulo 6

Solução. 6.3

library(mvtnorm)

m <- c(100,0,-10,-200)
rho <- matrix(c(1   ,-.86, .38,-.67,
                -.86,1   ,-.71, .25,
                .38 ,-.71,   1, .23,
                -.67, .25, .23,   1), byrow = TRUE, ncol = 4)
rho

##       [,1]  [,2]  [,3]  [,4]
## [1,]  1.00 -0.86  0.38 -0.67
## [2,] -0.86  1.00 -0.71  0.25
## [3,]  0.38 -0.71  1.00  0.23
## [4,] -0.67  0.25  0.23  1.00

rho2 <- rho
rho2[1,2] <- 0
rho2[1,3] <- 0
rho2[1,4] <- 0
rho2[2,3] <- 0
rho2[2,4] <- 0
rho2[3,4] <- 0

sigma <- rho
d <- c(3,16,18,61)
diag(sigma) <- d^2
sigma[1,2] <- sigma[2,1] <- sigma[1,2]*d[1]*d[2]
sigma[1,3] <- sigma[3,1] <- sigma[1,3]*d[1]*d[3]
sigma[1,4] <- sigma[4,1] <- sigma[1,4]*d[1]*d[4]
sigma[2,3] <- sigma[3,2] <- sigma[2,3]*d[2]*d[3]
sigma[2,4] <- sigma[4,2] <- sigma[2,4]*d[2]*d[4]
sigma[3,4] <- sigma[4,3] <- sigma[3,4]*d[3]*d[4]
sigma

##         [,1]    [,2]    [,3]    [,4]
## [1,]    9.00  -41.28   20.52 -122.61
## [2,]  -41.28  256.00 -204.48  244.00
## [3,]   20.52 -204.48  324.00  252.54
## [4,] -122.61  244.00  252.54 3721.00

# a
# https://stackoverflow.com/questions/47403003/why-is-the-pmvnorm-result-different-when-the-input-matrix-are-covariance-and-c

# b
?mvtnorm::pmvnorm

# c
# Pr(X1 < 100, X2 < 2, X3 < 0, X4 < -50)
lower <- c(-Inf, -Inf, -Inf, -Inf)
upper <- c(100, 2, 0, -50)
# GenzBretz
pmvnorm(lower, upper, mean = m, sigma=sigma)

## [1] 0.05095994
## attr(,"error")
## [1] 9.892229e-05
## attr(,"msg")
## [1] "Normal Completion"

pmvnorm(lower, upper = (upper-m)/sqrt(diag(sigma)), 
        mean = rep(0,nrow(sigma)), corr=rho)

## [1] 0.0509613
## attr(,"error")
## [1] 9.087318e-05
## attr(,"msg")
## [1] "Normal Completion"

# Miwa
pmvnorm(lower, upper, mean = m, sigma=sigma, algorithm = Miwa())

## [1] 0.05097291
## attr(,"error")
## [1] NA
## attr(,"msg")
## [1] "Normal Completion"

pmvnorm(lower, upper = (upper-m)/sqrt(diag(sigma)), 
        mean = rep(0,nrow(sigma)), corr=rho,
        algorithm = Miwa())

## [1] 0.05097291
## attr(,"error")
## [1] NA
## attr(,"msg")
## [1] "Normal Completion"

# TVPACK - somente para dimensões 2 ou 3

# d
# Pr(X_1 > 100, X_2 < 2, X_3 > 0, X_4 < -50)
lower <- c(100, -Inf, 0, -Inf)
upper <- c(Inf, 2, Inf, -50)

# GenzBretz
pmvnorm(lower, upper, mean = m, sigma=sigma)

## [1] 0.1965489
## attr(,"error")
## [1] 0.0001753258
## attr(,"msg")
## [1] "Normal Completion"

pmvnorm(lower = (lower-m)/sqrt(diag(sigma)), 
        upper = (upper-m)/sqrt(diag(sigma)), 
        mean = rep(0,nrow(sigma)), corr=rho)

## [1] 0.1966127
## attr(,"error")
## [1] 5.229961e-05
## attr(,"msg")
## [1] "Normal Completion"

# Miwa
pmvnorm(lower, upper, m, sigma=sigma, algorithm = Miwa())

## [1] 0.1965793
## attr(,"error")
## [1] NA
## attr(,"msg")
## [1] "Normal Completion"

pmvnorm(lower = (lower-m)/sqrt(diag(sigma)),
        upper = (upper-m)/sqrt(diag(sigma)), 
        mean = rep(0,nrow(sigma)), corr=rho,
        algorithm = Miwa())

## [1] 0.1965793
## attr(,"error")
## [1] NA
## attr(,"msg")
## [1] "Normal Completion"

# e
# Pr(70 < X_1 < 120, X_2 < 2, X_3 > 0, -125 < X_4 < -50)
lower <- c(70, -Inf, 0, -125)
upper <- c(120, 2, Inf, -50)

# GenzBretz
pmvnorm(lower, upper, mean = m, sigma=sigma)

## [1] 0.0313817
## attr(,"error")
## [1] 7.040005e-06
## attr(,"msg")
## [1] "Normal Completion"

pmvnorm(lower = (lower-m)/sqrt(diag(sigma)), 
        upper = (upper-m)/sqrt(diag(sigma)), 
        mean = rep(0,nrow(sigma)), corr=rho)

## [1] 0.03137679
## attr(,"error")
## [1] 4.720631e-06
## attr(,"msg")
## [1] "Normal Completion"

# Miwa
pmvnorm(lower, upper, m, sigma=sigma, algorithm = Miwa())

## [1] 0.03137884
## attr(,"error")
## [1] NA
## attr(,"msg")
## [1] "Normal Completion"

pmvnorm(lower = (lower-m)/sqrt(diag(sigma)),
        upper = (upper-m)/sqrt(diag(sigma)), 
        mean = rep(0,nrow(sigma)), corr=rho,
        algorithm = Miwa())

## [1] 0.03137884
## attr(,"error")
## [1] NA
## attr(,"msg")
## [1] "Normal Completion"

# f
set.seed(1) # para garantir reprodutibilidade
sim1 <- mvtnorm::rmvnorm(2000, mean = m, sigma = sigma)
head(sim1)

##           [,1]       [,2]       [,3]      [,4]
## [1,]  95.36375  15.545475 -19.402609 -104.2662
## [2,] 100.87836 -12.749509   6.736066 -156.9392
## [3,] 101.54520 -14.910552  18.333048 -173.0005
## [4,] 104.77871 -37.346878  23.510279 -205.3740
## [5,]  97.61704   9.066492  -1.546524 -157.5069
## [6,] 103.15600   1.183472 -21.040319 -319.5508

dim(sim1)

## [1] 2000    4

(ms <- colMeans(sim1)) # vetor de médias (amostrais)

## [1]   99.902383    0.478849  -10.374043 -199.559335

(covs <- cov(sim1))

##             [,1]       [,2]       [,3]      [,4]
## [1,]    9.319913  -42.12923   20.41825 -132.7420
## [2,]  -42.129228  258.67254 -207.36359  278.4131
## [3,]   20.418248 -207.36359  336.91002  261.8536
## [4,] -132.741956  278.41311  261.85363 3965.0054

(cors <- cor(sim1))

##            [,1]       [,2]       [,3]       [,4]
## [1,]  1.0000000 -0.8580299  0.3643809 -0.6905266
## [2,] -0.8580299  1.0000000 -0.7024257  0.2749113
## [3,]  0.3643809 -0.7024257  1.0000000  0.2265582
## [4,] -0.6905266  0.2749113  0.2265582  1.0000000

all.equal(cors, cov2cor(covs))

## [1] TRUE

# g
ICSNP::HotellingsT2(sim1, mu = m)

## 
##  Hotelling's one sample T2-test
## 
## data:  sim1
## T.2 = 1.0891, df1 = 4, df2 = 1996, p-value = 0.3602
## alternative hypothesis: true location is not equal to c(100,0,-10,-200)

# ?biotools::boxM

# h
# Pr(X1 < 100, X2 < 2, X3 < 0, X4 < -50)
lower <- c(-Inf, -Inf, -Inf, -Inf)
upper <- c(100, 2, 0, -50)
pmvnorm(lower, upper, mean = ms, sigma=covs)

## [1] 0.05032645
## attr(,"error")
## [1] 0.0001247989
## attr(,"msg")
## [1] "Normal Completion"

# Pr(X_1 > 100, X_2 < 2, X_3 > 0, X_4 < -50)
lower <- c(100, -Inf, 0, -Inf)
upper <- c(Inf, 2, Inf, -50)
pmvnorm(lower, upper, mean = ms, sigma=covs)

## [1] 0.1889455
## attr(,"error")
## [1] 3.466999e-05
## attr(,"msg")
## [1] "Normal Completion"

# Pr(70 < X_1 < 120, X_2 < 2, X_3 > 0, -125 < X_4 < -50)
lower <- c(70, -Inf, 0, -125)
upper <- c(120, 2, Inf, -50)
pmvnorm(lower, upper, mean = ms, sigma=covs)

## [1] 0.03170688
## attr(,"error")
## [1] 7.90073e-06
## attr(,"msg")
## [1] "Normal Completion"

# # i
# # Pr(X1 < 100, X2 < 2, X3 < 0, X4 < -50)
# lower <- c(-Inf, -Inf, -Inf, -Inf)
# upper <- c(100, 2, 0, -50)
# 
# normal <- pmvnorm(lower, upper, mean = m, sigma=sigma, seed = 1)
# t_inf_gl <- pmvt(lower, upper, ms, sigma=sigma, df = 0, seed = 1)