Capítulo 5

Solução. 5.3

library(mvtnorm)
library(matrixcalc)

m <- c(100, 0, -10, 200)
r12 <- -0.86 # correlação das variáveis 1 e 2
r13 <- 0.38
r14 <- -0.67
r23 <- -0.71
r24 <- 0.25
r34 <- 0.23
dp1 <- 3
dp2 <- 16
dp3 <- 18
dp4 <- 61
s <- matrix(c(dp1^2,       r12*dp1*dp2, r13*dp1*dp3, r14*dp1*dp4,
              r12*dp1*dp2, dp2^2,       r23*dp2*dp3, r24*dp2*dp4,
              r13*dp1*dp3, r23*dp2*dp3, dp3^2,       r34*dp3*dp4,
              r14*dp1*dp4, r24*dp2*dp4, r34*dp3*dp4, dp4^2),
            ncol = 4, nrow = 4, byrow = T)
colnames(s) <- paste0('X', 1:4)
rownames(s) <- paste0('X', 1:4)
s

##         X1      X2      X3      X4
## X1    9.00  -41.28   20.52 -122.61
## X2  -41.28  256.00 -204.48  244.00
## X3   20.52 -204.48  324.00  252.54
## X4 -122.61  244.00  252.54 3721.00

isSymmetric(s)

## [1] TRUE

matrixcalc::is.positive.definite(s)

## [1] TRUE

eigen(s)

## eigen() decomposition
## $values
## [1] 3758.6816969  500.7897519   50.2106611    0.3178901
## 
## $vectors
##             [,1]         [,2]        [,3]       [,4]
## [1,]  0.03287727 -0.087714383  0.12415761 0.98783104
## [2,] -0.06566280  0.651022239 -0.74055893 0.15307157
## [3,] -0.06904642 -0.753945967 -0.65306520 0.01743331
## [4,] -0.99490707  0.006458407  0.09830159 0.02133101

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

##           [,1]       [,2]       [,3]      [,4]
## [1,]  95.36375  15.545475 -19.402609 295.73379
## [2,] 100.87836 -12.749509   6.736066 243.06077
## [3,] 101.54520 -14.910552  18.333048 226.99947
## [4,] 104.77871 -37.346878  23.510279 194.62603
## [5,]  97.61704   9.066492  -1.546524 242.49308
## [6,] 103.15600   1.183472 -21.040319  80.44921

dim(sim1)

## [1] 2000    4

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

## [1]  99.902383   0.478849 -10.374043 200.440665

(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

# brindes
pairs(sim1)

# Pr(X1 < 0, X2 < 0, X3 < 0, X4 < 0)
lower <- c(-Inf, -Inf, -Inf, -Inf)
upper <- c(0, 0, 0, 0)
pmvnorm(lower, upper, ms, corr=cors)

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

pmvnorm(lower, upper, ms, sigma=covs)

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

# c
# Pr(X1 < 70, X2 < 1.96, X3 < 0, X4 < -50)
lower <- c(-Inf, -Inf, -Inf, -Inf)
upper <- c(70, 1.96, 0, -50)
pmvnorm(lower, upper, ms, corr=cors)

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

pmvnorm(lower, upper, ms, sigma=covs)

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

# The multivariate normal case is treated as a special case of pmvt with df=0 and univariate problems are passed to pnorm.
pmvt(lower, upper, ms, corr=cors, df = 0)

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

pmvt(lower, upper, ms, sigma=covs, df = 0)

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

# d
# Pr(X_1 > 70, X_2 < 1.96, X_3 > 0, X_4 < -50)
lower <- c(70, -Inf, 0, -Inf)
upper <- c(Inf, 1.96, Inf, -50)
pmvnorm(lower, upper, ms, corr=cors)

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

pmvnorm(lower, upper, ms, sigma=covs)

## [1] 2.085381e-06
## attr(,"error")
## [1] 7.209669e-10
## attr(,"msg")
## [1] "Normal Completion"

pmvt(lower, upper, ms, corr=cors, df = 0)

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

pmvt(lower, upper, ms, sigma=covs, df = 0)

## [1] 2.085731e-06
## attr(,"error")
## [1] 5.568954e-10
## attr(,"msg")
## [1] "Normal Completion"

# e
# Pr(70 < X_1 < 120, X_2 < 1.96, X_3 > 0, -125 < X_4 < -50)
lower <- c(70, -Inf, 0, -125)
upper <- c(120, 1.96, Inf, -50)
pmvnorm(lower, upper, ms, corr=cors)

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

pmvnorm(lower, upper, ms, sigma=covs)

## [1] 2.081597e-06
## attr(,"error")
## [1] 2.187142e-10
## attr(,"msg")
## [1] "Normal Completion"

pmvt(lower, upper, ms, corr=cors, df = 0)

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

pmvt(lower, upper, ms, sigma=covs, df = 0)

## [1] 2.08148e-06
## attr(,"error")
## [1] 2.265805e-10
## attr(,"msg")
## [1] "Normal Completion"