Resolução em R dos Exemplos do livro: Introdução a Modelos Lineares

João Gil de Luna
joaogill@yahoo.es
http://lattes.cnpq.br/9922128327962088

Ricardo Alves de Olinda
ricardo.estat@yahoo.com.br
http://lattes.cnpq.br/7767223263366578

Universidade Estadual da Paraíba - UEPB
Centro de Ciências e Tecnologia - CCT
Departamento de Estatística

Sumário

Carregando pacotes necessários

library(pracma)
library(MASS)
require(matrixcalc)

## Carregando pacotes exigidos: matrixcalc

require(Matrix)

## Carregando pacotes exigidos: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:pracma':
## 
##     expm, lu, tril, triu

require(ibd)

## Carregando pacotes exigidos: ibd

require(pracma)
require(matlib)

## Carregando pacotes exigidos: matlib

## 
## Attaching package: 'matlib'

## The following object is masked from 'package:matrixcalc':
## 
##     vec

## The following objects are masked from 'package:pracma':
## 
##     angle, inv

require(MASS)
require(tidyverse)

## Carregando pacotes exigidos: tidyverse

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ purrr::cross()  masks pracma::cross()
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ tidyr::pack()   masks Matrix::pack()
## ✖ dplyr::select() masks MASS::select()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

require(kableExtra)

## Carregando pacotes exigidos: kableExtra
## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows

require(DT)

## Carregando pacotes exigidos: DT

require(agricolae)

## Carregando pacotes exigidos: agricolae

Exemplos do Capítulo I

Tópicos de Matrizes

Exemplo 1.2.4

A <- matrix(c(2, 3, -1, 5, 1, 2), nrow = 2, byrow = TRUE)
A

##      [,1] [,2] [,3]
## [1,]    2    3   -1
## [2,]    5    1    2

At <- t(A) # calculando a transposta da matriz A
At

##      [,1] [,2]
## [1,]    2    5
## [2,]    3    1
## [3,]   -1    2

Exemplo 1.3.5

A <- matrix(c(2, -1, 0, 1, 2, -3), nrow = 2, byrow = TRUE)

B <- matrix(c(1, 1, 1, 1, 2, 3), nrow = 3, byrow = TRUE)

AB <- A %*% B
AB

##      [,1] [,2]
## [1,]    1    1
## [2,]   -3   -6

Exemplo 1.3.6

## A ⊕ B =

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    2    3    0    0
## [2,]    0    0    0    6    7
## [3,]    0    0    0    4   -1

Exemplo 1.3.7

## A ⊗ v =

##      [,1] [,2]
## [1,]    1    0
## [2,]    2    0
## [3,]    3    0
## [4,]    0    1
## [5,]    0    2
## [6,]    0    3

## 
## v ⊗ A =

##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
## [3,]    2    0
## [4,]    0    2
## [5,]    3    0
## [6,]    0    3

Exemplo 1.3.8

## A^2 =

##      [,1] [,2]
## [1,]    7   10
## [2,]   15   22

## 
## A^3 =

##      [,1] [,2]
## [1,]   37   54
## [2,]   81  118

Exemplo 1.6.2

## Forma escalonada canônica de A:

##      [,1] [,2] [,3]
## [1,]    1    0    1
## [2,]    0    1   -1
## [3,]    0    0    0

Exemplo 1.7.2

## [1] "Matriz aumentada [A | I]:"

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    4    2    2    1    0    0
## [2,]    2    2    0    0    1    0
## [3,]    2    0    2    0    0    1

## 
## Após eliminar coluna 1:

##      [,1] [,2] [,3]  [,4] [,5] [,6]
## [1,]    1  0.5  0.5  0.25    0    0
## [2,]    0  1.0 -1.0 -0.50    1    0
## [3,]    0 -1.0  1.0 -0.50    0    1

## 
## Após eliminar coluna 2:

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    0    1  0.5 -0.5    0
## [2,]    0    1   -1 -0.5  1.0    0
## [3,]    0    0    0 -1.0  1.0    1

## 
## Matriz final (escalonada):

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    0    1  0.5 -0.5    0
## [2,]    0    1   -1 -0.5  1.0    0
## [3,]    0    0    0  1.0 -1.0   -1

Exemplo 1.8.1

## [1] "Matriz A:"

##      [,1] [,2] [,3]
## [1,]    1    2    1
## [2,]    2    1    0
## [3,]    1    1    1

## [1] "Matriz Inversa de A:"

##      [,1] [,2] [,3]
## [1,] -0.5  0.5  0.5
## [2,]  1.0  0.0 -1.0
## [3,] -0.5 -0.5  1.5

## [1] "Verificação A %*% A_inv:"

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

Exemplo 1.9.1

## Matriz aumentada [A | I]:

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    2    1    1    0    0
## [2,]    2    3    1    0    1    0
## [3,]    1    1    1    0    0    1

## 
## Após zerar abaixo da diagonal:

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    2    1    1    0    0
## [2,]    0   -1   -1   -2    1    0
## [3,]    0   -1    0   -1    0    1

## 
## Após L3 <- L3 + L2:

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    2    1    1    0    0
## [2,]    0   -1   -1   -2    1    0
## [3,]    0   -2   -1   -3    1    1

## 
## Matriz [D | C] final:

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    2    1    1    0    0
## [2,]    0    1    1    2   -1    0
## [3,]    0   -2   -1   -3    1    1

## 
## Matriz D (diagonal):

##      [,1] [,2] [,3]
## [1,]    1    2    1
## [2,]    0    1    1
## [3,]    0   -2   -1

## 
## Matriz C (transformação por congruência):

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    2   -1    0
## [3,]   -3    1    1

## 
## Verificação D = C A C':

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0   -1    2
## [3,]    0    2   -3

Exemplo 1.9.3

## Autovalores (D):

##      [,1] [,2] [,3]
## [1,]    4    0    0
## [2,]    0    1    0
## [3,]    0    0    1

## Autovetores ajustados (C):

##        [,1]    [,2]    [,3]
## [1,] 0.5774  0.0000  0.8165
## [2,] 0.5774 -0.7071 -0.4082
## [3,] 0.5774  0.7071 -0.4082

## C^-1:

##        [,1]    [,2]    [,3]
## [1,] 0.5774  0.5774  0.5774
## [2,] 0.0000 -0.7071  0.7071
## [3,] 0.8165 -0.4082 -0.4082

## Raiz quadrada de D:

##      [,1] [,2] [,3]
## [1,]    2    0    0
## [2,]    0    1    0
## [3,]    0    0    1

Exemplo 1.10.1

## [1] 3 1

Exemplo 1.11.2

##    u1 u2
## u1  3  0
## u2  0  1

Exemplo 1.12.1

##      [,1] [,2]          [,3]
## [1,]    2    0  3.441276e-08
## [2,]    1   -1 -3.441276e-08
## [3,]    1    1 -3.441276e-08

##      [,1]         [,2]         [,3]
## [1,]    4 2.000000e+00 2.000000e+00
## [2,]    2 2.000000e+00 1.184238e-15
## [3,]    2 1.184238e-15 2.000000e+00

Exemplos do Capítulo II

Soluções de Equações Lineares

Exemplo 2.2.1

O sistema de equações lineares \(Ax = g\) é dado por:

\[ \begin{cases} 4x_1 + 2x_2 + 2x_3 = 14 \\ 2x_1 + 2x_2 = 6 \\ 2x_1 + 2x_3 = 8 \end{cases} \]

Ou, na forma matricial:

\[ A = \begin{pmatrix} 4 & 2 & 2 \\ 2 & 2 & 0 \\ 2 & 0 & 2 \end{pmatrix}, \quad x = \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix}, \quad g = \begin{pmatrix} 14 \\ 6 \\ 8 \end{pmatrix} \]

A <- matrix(c(4, 2, 2,
              2, 2, 0,
              2, 0, 2), nrow = 3, byrow = TRUE)

g <- c(14, 6, 8)

A título de ilustração considere três vetores soluções.

Vetor solução \(x^{(1)} = \begin{pmatrix} 0 \\ 3 \\ 4 \end{pmatrix}\)

x1 <- c(0, 3, 4)
A %*% x1

##      [,1]
## [1,]   14
## [2,]    6
## [3,]    8

Vetor solução \(x^{(2)} = \begin{pmatrix} 3 \\ 0 \\ 1 \end{pmatrix}\)

x2 <- c(3, 0, 1)
A %*% x2

##      [,1]
## [1,]   14
## [2,]    6
## [3,]    8

Vetor solução \(x^{(3)} = \begin{pmatrix} 4 \\ -1 \\ 0 \end{pmatrix}\)

x3 <- c(4, -1, 0)
A %*% x3

##      [,1]
## [1,]   14
## [2,]    6
## [3,]    8

Todos os vetores \(x^{(j)}\) satisfazem \(Ax = g\), evidenciando que o sistema é consistente e indeterminado.

Exemplo 2.2.2

Considerando o sistema de equações lineares \(Ax = g\) do exemplo anterior:

\[ A = \begin{pmatrix} 4 & 2 & 2 \\ 2 & 2 & 0 \\ 2 & 0 & 2 \end{pmatrix}, \quad g = \begin{pmatrix} 14 \\ 6 \\ 8 \end{pmatrix} \]

A <- matrix(c(4, 2, 2,
              2, 2, 0,
              2, 0, 2), nrow = 3, byrow = TRUE)

g <- matrix(c(14, 6, 8), ncol = 1)

A_aumentada <- cbind(A, g)
A_aumentada

##      [,1] [,2] [,3] [,4]
## [1,]    4    2    2   14
## [2,]    2    2    0    6
## [3,]    2    0    2    8

rref_Ag <- rref(A_aumentada)
rref_Ag

##      [,1] [,2] [,3] [,4]
## [1,]    1    0    1    4
## [2,]    0    1   -1   -1
## [3,]    0    0    0    0

rank_A <- qr(A)$rank
rank_Ag <- qr(A_aumentada)$rank

rank_A

## [1] 2

rank_Ag

## [1] 2

A análise da matriz aumentada confirma que: r(A) = r(A:g) = 2

Portanto, o sistema é consistente, como já verificado no Exemplo 2.2.1.

Exemplo 2.2.3

Suponhamos agora o sistema \[ \begin{pmatrix} 4 & 2 & 2 \\ 2 & 2 & 0 \\ 2 & 0 & 2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix} \]

A <- matrix(c(4, 2, 2,
              2, 2, 0,
              2, 0, 2), nrow = 3, byrow = TRUE)

g <- matrix(c(1, 1, 1), ncol = 1)

A_aumentada <- cbind(A, g)
A_aumentada

##      [,1] [,2] [,3] [,4]
## [1,]    4    2    2    1
## [2,]    2    2    0    1
## [3,]    2    0    2    1

rref_Ag <- rref(A_aumentada)
rref_Ag

##      [,1] [,2] [,3] [,4]
## [1,]    1    0    1    0
## [2,]    0    1   -1    0
## [3,]    0    0    0    1

rank_A <- qr(A)$rank
rank_Ag <- qr(A_aumentada)$rank

rank_A

## [1] 2

rank_Ag

## [1] 3

\(r(A) = 2 \ne r(A:g) = 3\) logo o sistema é inconsistente.

Exemplo 2.2.4

Seja o sistema \(Ax = g\) caracterizado por:

\[ \begin{pmatrix} 2 & 1\\ 1 & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 7 \\ 5 \end{pmatrix} \]

Então, \(x^0\) será:

A <- matrix(c(2, 1,
              1, 1), nrow = 2, byrow = TRUE)

g <- c(7, 5)

x_sol <- solve(A, g)
x_sol

## [1] 2 3

Note que nesse caso a solução é única. Isto é, existe apenas uma combinação linear das colunas de A que reproduz g. Ou seja,

c1 <- A[,1]
c2 <- A[,2]

g_verif <- 2 * c1 + 3 * c2
g_verif

## [1] 7 5

Nesse caso o Teorema 2.2.2 pode fornecer diretamente a solução.

A_aumentada <- cbind(A, g)
rref_Ag <- rref(A_aumentada)
rref_Ag

##          g
## [1,] 1 0 2
## [2,] 0 1 3

Exemplo 2.2.5

Seja o sistema \(Ax = g\) caracterizado por:

\[ \begin{pmatrix} 1 & 1 \\ 1 & -1 \\ -2 & 0 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 3 \\ 1 \\ 2 \end{pmatrix} \]

Usando o algoritmo de Searle, temos

\[ M = \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} e \quad M^{-1} = \begin{pmatrix} \frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & -\frac{1}{2} \end{pmatrix} \]

Portanto,

\[ A^- = \begin{pmatrix} \frac{1}{2} & \frac{1}{2} & 0 \\ \frac{1}{2} & -\frac{1}{2} & 0 \end{pmatrix} \]

É fácil verificar que \(AA^-g \ne g\)

A <- matrix(c(1, 1,
              1, -1,
              -2, 0), nrow = 3, byrow = TRUE)

g <- matrix(c(3, 1, 2), ncol = 1)

A_inv <- matrix(c(1/2, 1/2, 0,
                  1/2, -1/2, 0), nrow = 2, byrow = TRUE)

A %*% A_inv %*% g

##      [,1]
## [1,]    3
## [2,]    1
## [3,]   -4

Portanto o sistema é inconsistente.

Exemplo 2.2.6

Considere o sistema \(Ax = g\) caracterizado por:

\[ \begin{pmatrix} 1 & 1 \\ 1 & -1 \\ -2 & 0 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 3 \\ 1 \\ -4 \end{pmatrix} \]

Usando o fato de que

\[ A^- = \begin{pmatrix} \frac{1}{2} & \frac{1}{2} & 0 \\ \frac{1}{2} & -\frac{1}{2} & 0 \end{pmatrix}, \quad A^\ell = \frac{1}{6} \begin{pmatrix} 1 & 1 & -2 \\ 3 & -3 & 0 \end{pmatrix} = A^+ \]

Tem-se que

\[ AA^-g = AA^\ell g = AA^+g = \begin{pmatrix} 3 \\ 1 \\ -4 \end{pmatrix} = g \]

e o sistema é consistente.

A <- matrix(c(1, 1,
              1, -1,
              -2, 0), nrow = 3, byrow = TRUE)

g <- matrix(c(3, 1, -4), ncol = 1)

A_inv <- matrix(c(1/2, 1/2, 0,
                  1/2, -1/2, 0), nrow = 2, byrow = TRUE)

A_ell <- (1/6) * matrix(c(1, 1, -2,
                          3, -3, 0), nrow = 2, byrow = TRUE)

A %*% A_inv %*% g

##      [,1]
## [1,]    3
## [2,]    1
## [3,]   -4

A %*% A_ell %*% g

##      [,1]
## [1,]    3
## [2,]    1
## [3,]   -4

Exemplo 2.2.7

Seja \(Ax = g\) consistente, dado por

\[ \begin{pmatrix} 4 & 2 & 2 \\ 2 & 2 & 0 \\ 2 & 0 & 2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 14 \\ 6 \\ 8 \end{pmatrix} \]

Tomando três inversas condicionais de A, temos três soluções distintas para o sistema.

g <- matrix(c(14, 6, 8), ncol = 1)

# Inversas condicionais
A1_inv <- matrix(c(0, 0, 0,
                   0, 1/2, 0,
                   0, 0, 1/2), nrow = 3, byrow = TRUE)

A2_inv <- matrix(c(1/2, -1/2, 0,
                   -1/2, 1, 0,
                   0, 0, 0), nrow = 3, byrow = TRUE)

A3_inv <- matrix(c(1/2, 0, -1/2,
                   0, 0, 0,
                   -1/2, 0, 1), nrow = 3, byrow = TRUE)

# Soluções
x1 <- A1_inv %*% g
x2 <- A2_inv %*% g
x3 <- A3_inv %*% g

x1

##      [,1]
## [1,]    0
## [2,]    3
## [3,]    4

x2

##      [,1]
## [1,]    4
## [2,]   -1
## [3,]    0

x3

##      [,1]
## [1,]    3
## [2,]    0
## [3,]    1

Tomemos agora o vetor:

\[ x^0 = A_1^- g + (I - A_1^- A) h = \begin{pmatrix} 0 \\ 3 \\ 4 \end{pmatrix} + \begin{pmatrix} 1 & 0 & 0 \\ -1 & 0 & 0 \\ -1 & 0 & 0 \end{pmatrix} \begin{pmatrix} h_1 \\ h_2 \\ h_3 \end{pmatrix} \]

Portanto,

\[ x^0 = \begin{pmatrix} h_1 \\ 3 - h_1 \\ 4 - h_1 \end{pmatrix} \quad \text{é solução para todo } h_1 \in \mathbb{R}. \]

Exemplo 2.4.1

A solução do sistema inconsistente \(X\theta = y\), do exemplo anterior, dada por:

X <- matrix(c(1, 1, 0,
              1, 1, 0,
              1, 0, 1,
              1, 0, 1), nrow = 4, byrow = TRUE)

y <- c(2, 3, 5, 4)

# Pseudoinversa de X
X_pinv <- ginv(X)

# Melhor solução aproximada
theta <- X_pinv %*% y
theta

##           [,1]
## [1,] 2.3333333
## [2,] 0.1666667
## [3,] 2.1666667

é a melhor solução aproximada.

Exemplo 2.5.3

Para o exemplo em questão, temos

y <- c(2, 3, 5, 4)

y_hat <- 0.5 * c(5, 5, 9, 9)

e_hat <- y - y_hat
e_hat

## [1] -0.5  0.5  0.5 -0.5

Exemplo 2.5.4

Para o exemplo em questão, é imediato verificar que

\[ X'ê = 0. \]

Ou seja,

X <- matrix(c(1, 1, 1, 1,
               1, 1, 0, 0,
               0, 0, 1, 1), nrow = 3, byrow = TRUE)

X %*% e_hat

##      [,1]
## [1,]    0
## [2,]    0
## [3,]    0

Exemplos do Capítulo III

Classificação das Formas Quadráticas

Exemplo 3.3.1

Pode-se verificar a classificação da forma quadrática atráves da classificação da matriz núcleo. Ou seja, uma forma quadrática tem a mesma classificação de sua matriz matriz núcleo. Para isso basta diagonalizar a matriz núcleo atráves de operações de congruência.

Por exemplo: \(A = \begin{pmatrix} 1 & 3 \\ 3 & 1 \end{pmatrix}\). Então,

A <- matrix(c(1,3,3,1), 2, 2, byrow = TRUE)

# Operação de linha (L2 <- L2 - 3*L1)
A[2, ] <- A[2, ] - 3*A[1, ]

# Mesma operação na coluna correspondente (C2 <- C2 - 3*C1)
A[ ,2] <- A[ ,2] - 3*A[ ,1]

A

##      [,1] [,2]
## [1,]    1    0
## [2,]    0   -8

Visto que a diagonal possui um valor positivo e um valor negativo conclui-se que a matriz núcleo de A é não definida, logo A é não definida.

Exemplo 3.3.2

Seja
\[ Q(x) = 2x_1^2 + 2x_1x_2 + 3x_2^2 \]

Portanto, \(A = \begin{pmatrix} 2 & 1 \\ 1 & 3 \end{pmatrix}\) e sua diagonal congruente \(D_1 = \begin{pmatrix} 2 & 0 \\ 0 & \frac{5}{2} \end{pmatrix}\).

Portanto, \(Q(x)\) é positiva definida.

Tomando como exemplo, o vetor \(x = \begin{pmatrix} 10 \\ 20 \end{pmatrix}\). Então:

A <- matrix(c(2, 1,
              1, 3), 
            nrow = 2, byrow = TRUE)

x <- c(10, 20)

Qx <- t(x) %*% A %*% x
Qx

##      [,1]
## [1,] 1800

Seja, agora, a transformação não singular \(x = Cy\), onde:

C <- matrix(c(2, 3,
              4, 5),
            nrow = 2, byrow = TRUE)
C

##      [,1] [,2]
## [1,]    2    3
## [2,]    4    5

e a inversa de C

C_inv <- solve(C)
C_inv

##      [,1] [,2]
## [1,] -2.5  1.5
## [2,]  2.0 -1.0

Temos \(y = C^{-1}x\) e \(M = C{´}AC\)

y <- C_inv %*% x
y

##      [,1]
## [1,]    5
## [2,]    0

M <- t(C) %*% A %*% C
M

##      [,1] [,2]
## [1,]   72   94
## [2,]   94  123

Como \(Q(y) = y^{´}My\) então:

Qy <- t(y) %*% M %*% y
Qy

##      [,1]
## [1,] 1800

Exemplo 3.3.3

Seja

\[ Q(\mathbf{x}) = \mathbf{x}^\top A \mathbf{x}, \quad A = \begin{pmatrix} 2 & -1 & -1 \\ -1 & 2 & -1 \\ -1 & -1 & 2 \end{pmatrix} \]

e seja a transformação \(x = Uy\), em que \(U^{´} = \begin{pmatrix} \frac{-2}{\sqrt 6} & \frac{1}{\sqrt 6} & \frac{1}{\sqrt 6} \\ 0 & \frac{1}{\sqrt 2} & \frac{-1}{\sqrt 3} \\ \frac{1}{\sqrt 3} & \frac{1}{\sqrt 3} & \frac{1}{\sqrt 3} \end{pmatrix}\). Então:

A <- matrix(c(2, -1, -1,
              -1, 2, -1,
              -1, -1, 2), nrow = 3, byrow = TRUE)

decomp <- eigen(A)
vetores <- decomp$vectors
U <- vetores

Lambda <- t(U) %*% A %*% U
zapsmall(Lambda)

##      [,1] [,2] [,3]
## [1,]    3    0    0
## [2,]    0    3    0
## [3,]    0    0    0

Assim, \(\lambda_1 = 3, \lambda_2 = 3\) e \(\lambda_3 = 0\). Portanto, \(Q(x)\) é semi positiva definida.

Tomando, como exemplo, o vetor \(x = \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}\). Então:

Q <- function(x) t(x) %*% A %*% x
x <- c(1, 1, 1)
Q(x)

##      [,1]
## [1,]    0

y <- t(U) %*% x
zapsmall(y)

##          [,1]
## [1,] 0.000000
## [2,] 0.000000
## [3,] 1.732051

Qy <- t(y) %*% Lambda %*% y
Qy[abs(Qy) < 1e-12] <- 0
Qy

##      [,1]
## [1,]    0

Naturalmente, dada a similaridade, \(Q(x)\) e \(Q(y)\) tem a mesma classificação.

Exemplo 3.3.4

Seja a forma quadrática \(Q(x) = x^{´}Ax\), onde

A <- matrix(c(3, 1, 1,
              1, 3, 1,
              1, 1, 3), nrow = 3, byrow = TRUE)

decomp <- eigen(A)
valores <- decomp$values

Lambda <- diag(valores)
Lambda

##      [,1] [,2] [,3]
## [1,]    5    0    0
## [2,]    0    2    0
## [3,]    0    0    2

Assim, \(\lambda_1 = 3, \lambda_2 = 3\) e \(\lambda_3 = 0\), portanto, \(A\) e \(Q(x)\) é positiva definida.

Seja, como ilustração \(x = \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}\)

Q <- function(x) t(x) %*% A %*% x

x <- c(1, 1, 1)

Q(x)

##      [,1]
## [1,]   15

Por outro lado, para \(y = B^{´}x\) onde \(B = C^{-1}D^{\frac{1}{2}}\)

vetores <- decomp$vectors 
B <- vetores
y <- t(B) %*% x

Qy <- t(y) %*% Lambda %*% y
Qy

##      [,1]
## [1,]   15

Exemplos do Capítulo IV

Introdução aos Modelos Lineares

Exemplo 4.3.1

Associando o modelo \(\mathbf{y}=\mathbf{X}\boldsymbol{\theta}+\mathbf{e}\) caracterizado por \(y_{ij}=\mu+\tau_i+e_{ij}\) às observações obtidas no experimento de Exemplo 4.3.1, podemos escrever:

\[ \begin{pmatrix} y_{11} \\ y_{12} \\ y_{13} \\ y_{14} \\ y_{21} \\ y_{22} \\ y_{23} \\ y_{31} \\ y_{32} \\ y_{33} \\ \end{pmatrix} = \begin{pmatrix} 5,0 \\ 4,0 \\ 3,0 \\ 4,0 \\ 6,0 \\ 7,0 \\ 8,0 \\ 9,0 \\ 8,0 \\ 10,0 \\ \end{pmatrix}= \begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \\ \end{pmatrix} \begin{pmatrix} \mu \\ \tau_1 \\ \tau_2 \\ \tau_3 \\ \end{pmatrix}+ \begin{pmatrix} e_{11} \\ e_{12} \\ e_{13} \\ e_{14} \\ e_{21} \\ e_{22} \\ e_{23} \\ e_{31} \\ e_{32} \\ e_{33} \\ \end{pmatrix} \]

Assim, para os dados do Exemplo 4.3.1, temos o seguinte S.E.N. \(\mathbf{X^{'}}\mathbf{X}\boldsymbol{\theta^{o}}=\mathbf{X^{'}}\mathbf{y}\)

• Implementação no R

X=matrix(c(1, 1, 0, 0,1, 1, 0, 0,
1, 1, 0, 0,1, 1, 0, 0,
1, 0, 1, 0,1, 0, 1, 0,
1, 0, 1, 0,1, 0, 0, 1,
1, 0, 0, 1,1, 0, 0, 1),byrow=T,nrow=10)
y=matrix(c(5,4,3,4,6,
           7,8,9,8,10),byrow=T,nrow=10)
showEqn(X,y,paste0("beta", 1:ncol(X)))

## 1*beta1 + 1*beta2 + 0*beta3 + 0*beta4  =   5 
## 1*beta1 + 1*beta2 + 0*beta3 + 0*beta4  =   4 
## 1*beta1 + 1*beta2 + 0*beta3 + 0*beta4  =   3 
## 1*beta1 + 1*beta2 + 0*beta3 + 0*beta4  =   4 
## 1*beta1 + 0*beta2 + 1*beta3 + 0*beta4  =   6 
## 1*beta1 + 0*beta2 + 1*beta3 + 0*beta4  =   7 
## 1*beta1 + 0*beta2 + 1*beta3 + 0*beta4  =   8 
## 1*beta1 + 0*beta2 + 0*beta3 + 1*beta4  =   9 
## 1*beta1 + 0*beta2 + 0*beta3 + 1*beta4  =   8 
## 1*beta1 + 0*beta2 + 0*beta3 + 1*beta4  =  10

• Uma solução do sistema de equações:

theta1=c("mu","tau1","tau2","tau3")
theta1

## [1] "mu"   "tau1" "tau2" "tau3"

theta_bola=matrix(c(0,4,7,9),byrow=T,nrow=4)
yhat=X%*%theta_bola
yhat

##       [,1]
##  [1,]    4
##  [2,]    4
##  [3,]    4
##  [4,]    4
##  [5,]    7
##  [6,]    7
##  [7,]    7
##  [8,]    9
##  [9,]    9
## [10,]    9

• Para os dados do Exemplo 4.3.1, se considerarmos a inversa condicional mais simples de \(\mathbf{X^{'}}\mathbf{X}\), teremos

InvXtX1=matrix(c(0, 0, 0,0,
                0,1/4,0,0,
                0,0,1/3,0,
                0,0,0, 1/3),byrow=T,nrow=4)
InvXtX1

##      [,1] [,2]      [,3]      [,4]
## [1,]    0 0.00 0.0000000 0.0000000
## [2,]    0 0.25 0.0000000 0.0000000
## [3,]    0 0.00 0.3333333 0.0000000
## [4,]    0 0.00 0.0000000 0.3333333

• Análise da Variância

Input =("
obs;tratamento;peso
1;T1;5
2;T1;4
3;T1;3
4;T1;4
5;T2;6
6;T2;7
7;T2;8
8;T3;9
9;T3;8
10;T3;10
")

dados=read.table(textConnection(Input),header=TRUE,sep=";")
ordem_tratamentos=c("T1","T2","T3")
dados$tratamento=factor(dados$tratamento,levels=ordem_tratamentos)
dados = dados %>%
  arrange(tratamento)

datatable(dados,options = list(
            scrollX = TRUE,
            scrollY = "400px",
            pageLength = 10,
            paging = TRUE
          )
)

• Visualização gráfica

dados %>% 
  ggplot(aes(x=tratamento,y=peso,fill=tratamento)) +
         geom_boxplot()

• Apresentação de forma matricial

anova=aov(peso ~ tratamento , data=dados)
X=as.matrix(model.matrix(anova))
y=matrix(dados$peso,ncol=1)
showEqn(X,y,c("mu","tau1","tau2","tau3"))

## 1*mu + 0*tau1 + 0*tau2  =   5 
## 1*mu + 0*tau1 + 0*tau2  =   4 
## 1*mu + 0*tau1 + 0*tau2  =   3 
## 1*mu + 0*tau1 + 0*tau2  =   4 
## 1*mu + 1*tau1 + 0*tau2  =   6 
## 1*mu + 1*tau1 + 0*tau2  =   7 
## 1*mu + 1*tau1 + 0*tau2  =   8 
## 1*mu + 0*tau1 + 1*tau2  =   9 
## 1*mu + 0*tau1 + 1*tau2  =   8 
## 1*mu + 0*tau1 + 1*tau2  =  10

• Sistema de Equações Normais:

showEqn(t(X)%*%X,t(X)%*%y,c("mu","tau1","tau2","tau3"))

## 10*mu + 3*tau1 + 3*tau2  =  64 
##  3*mu + 3*tau1 + 0*tau2  =  21 
##  3*mu + 0*tau1 + 3*tau2  =  27

Para o exemplo que estamos considerando, temos

SQTot=t(y)%*%y
SQTot

##      [,1]
## [1,]  460

• Soma de Quadrado dos Parâmetros

P=X%*%solve(t(X)%*%X)%*%t(X)
SQPar=t(y)%*%P%*%y
SQPar

##      [,1]
## [1,]  454

• Soma de Quadrado dos Resíduos

SQRes=SQTot-SQPar
SQRes

##      [,1]
## [1,]    6

• Soma de Quadrado de Tratamentos

X1=X[,1]
P1=X1%*%solve(t(X1)%*%X1)%*%t(X1)
SQTrat=t(y)%*%(P-P1)%*%y
SQTrat

##      [,1]
## [1,] 44.4

que confere com a Soma de Quadrados de Tratamentos da ANOVA

anova1= aov(peso~tratamento, data=dados)
summary(anova1)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## tratamento   2   44.4  22.200    25.9 0.000582 ***
## Residuals    7    6.0   0.857                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Definir contrastes ortogonais (exemplo)
contrastes <- matrix(c(0, -1, 1,
                      -2,  1, 1),
                     nrow = 2,
                     byrow = TRUE)
colnames(contrastes) <- levels(dados$tratamento)

contrasts(dados$tratamento)

##    T2 T3
## T1  0  0
## T2  1  0
## T3  0  1

mat.c <- matrix(c(0,-1,1, -2,1,1), nc=2)
mat.c

##      [,1] [,2]
## [1,]    0   -2
## [2,]   -1    1
## [3,]    1    1

contrasts(dados$tratamento) <- mat.c
summary(anova1, split=list("tratamento"=list("C.orto.1"=1, "C.orto.2"=2)))

##                        Df Sum Sq Mean Sq F value   Pr(>F)    
## tratamento              2  44.40   22.20    25.9 0.000582 ***
##   tratamento: C.orto.1  1   1.54    1.54     1.8 0.221601    
##   tratamento: C.orto.2  1  42.86   42.86    50.0 0.000199 ***
## Residuals               7   6.00    0.86                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estimação por intervalo

Exemplo 4.7.1

Exemplo 4.7.1 Encontre intervalos com \(99\%\) de confiança para as funções paramétricas \(\boldsymbol{\lambda}^{'}_1\boldsymbol{\theta}=-\tau_2+\tau_3\) e \(\boldsymbol{\lambda}^{'}_2\boldsymbol{\theta}=-2\tau_1+\tau_2+\tau_3\)

Com base nos dados do Exemplo 4.6.2, temos que

lambda1=c(0,0,-1,1)
lambda1_hat=lambda1%*%theta_bola
lambda1_hat

##      [,1]
## [1,]    2

e

lambda2=c(0,-2,1,1)
lambda2_hat=lambda2%*%theta_bola
lambda2_hat

##      [,1]
## [1,]    8

• Para um nível de confiança de 99% e 7 graus de liberdade (teste bilateral)

alpha <- 0.01
graus_liberdade <- 7
valor_t <- qt(1 - alpha/2, df = graus_liberdade)
print(valor_t)

## [1] 3.499483

• Por outro lado

estima1=lambda1%*%InvXtX1%*%lambda1

e

estima2=lambda2%*%InvXtX1%*%lambda2

anova2<-anova(anova1)
s2=anova2$Mean[2]

limite_inferior1 <- lambda1_hat - valor_t*sqrt(estima1*s2)
limite_superior1 <- lambda1_hat + valor_t*sqrt(estima1*s2)

resultado1 <- data.frame(limite_inferior = limite_inferior1,
                         limite_superior = limite_superior1)
print(resultado1)

##   limite_inferior limite_superior
## 1      -0.6453607        4.645361

limite_inferior2 <- lambda2_hat - valor_t*sqrt(estima2*s2)
limite_superior2 <- lambda2_hat + valor_t*sqrt(estima2*s2)

resultado2 <- data.frame(limite_inferior = limite_inferior2,
                         limite_superior = limite_superior2)
print(resultado2)

##   limite_inferior limite_superior
## 1        3.817317        12.18268

Exemplos do Capítulo V

Aplicações em Modelos Lineares

Exemplo 5.1 Considerando-se os dados da Tabela 5.3, os quais se referem à pureza do oxigênio \(y\) e a porcentagem de hidrocarbonetos presentes na unidade principal de destilação \(x\) observados num processo químico de destilação.

hidrocarboneto=c(0.99,1.02,1.15,1.29,1.46,1.36,0.87,1.23,1.55,1.40, 
1.19,1.15,0.98,1.01,1.11,1.20,1.26,1.32,1.43,0.95)

oxigênio=c(90.01,89.05,91.43,93.74,96.73,94.45,87.59,91.77,99.42,93.65,
  93.54,92.52,90.56,89.54,89.85,90.39,93.25,93.41,94.98,87.33)

# Combina os vetores em um data frame
pureza <- data.frame(hidrocarboneto, oxigênio)
pureza

##    hidrocarboneto oxigênio
## 1            0.99    90.01
## 2            1.02    89.05
## 3            1.15    91.43
## 4            1.29    93.74
## 5            1.46    96.73
## 6            1.36    94.45
## 7            0.87    87.59
## 8            1.23    91.77
## 9            1.55    99.42
## 10           1.40    93.65
## 11           1.19    93.54
## 12           1.15    92.52
## 13           0.98    90.56
## 14           1.01    89.54
## 15           1.11    89.85
## 16           1.20    90.39
## 17           1.26    93.25
## 18           1.32    93.41
## 19           1.43    94.98
## 20           0.95    87.33

attach(pureza)

## The following objects are masked _by_ .GlobalEnv:
## 
##     hidrocarboneto, oxigênio

• Algumas estatísticas descritivas

summary(hidrocarboneto)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.870   1.018   1.195   1.196   1.330   1.550

summary(oxigênio)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   87.33   89.97   92.14   92.16   93.67   99.42

cor(hidrocarboneto,oxigênio)

## [1] 0.9367154

• Ajuste do modelo de Regressão

X=matrix(c(1,0.99, 1,1.02,
           1,1.15, 1,1.29,
           1,1.46, 1,1.36,
           1,0.87, 1,1.23, 
           1,1.55, 1,1.40, 
           1,1.19, 1,1.15,
           1,0.98, 1,1.01,
           1,1.11, 1,1.20,
           1,1.26, 1,1.32,
           1,1.43, 1,0.95),byrow=T,nrow=20)
y=matrix(c(90.01,89.05,91.43,93.74,96.73,94.45,87.59,91.77,99.42,93.65,
  93.54,92.52,90.56,89.54,89.85,90.39,93.25,93.41,94.98,87.33),byrow=T,nrow=20)
y

##        [,1]
##  [1,] 90.01
##  [2,] 89.05
##  [3,] 91.43
##  [4,] 93.74
##  [5,] 96.73
##  [6,] 94.45
##  [7,] 87.59
##  [8,] 91.77
##  [9,] 99.42
## [10,] 93.65
## [11,] 93.54
## [12,] 92.52
## [13,] 90.56
## [14,] 89.54
## [15,] 89.85
## [16,] 90.39
## [17,] 93.25
## [18,] 93.41
## [19,] 94.98
## [20,] 87.33

• Estimação do vetor de parâmetros \(\boldsymbol{\beta}\)

XtX=t(X)%*%X
Inv_XtX=solve(XtX)
beta_hat=Inv_XtX%*%t(X)%*%y
beta_hat

##          [,1]
## [1,] 74.28331
## [2,] 14.94748

• O modelo de regressão linear simples ajustado é:

\(\hat{y}=74.2833+14.9475x\)

• Análise da Variância

fit=lm(oxigênio~hidrocarboneto)
anova(fit)

## Analysis of Variance Table
## 
## Response: oxigênio
##                Df Sum Sq Mean Sq F value    Pr(>F)    
## hidrocarboneto  1 152.13 152.127  128.86 1.227e-09 ***
## Residuals      18  21.25   1.181                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Variação Total

pureza$dif <- pureza$oxigênio - mean(pureza$oxigênio)

• O intercepto e a inclinação da equação da reta ajustada:

(cf <- coef(fit))

##    (Intercept) hidrocarboneto 
##       74.28331       14.94748

• Valores de \(Y\) previstos pela equação da reta para cada valor de \(x\):

pureza$Y.pred <- predict(fit)

• A diferença entre os valores de \(Y\) e os previstos, que são os resíduos da regressão:

pureza$residuo <- pureza$oxigênio - pureza$Y.pred
pureza

##    hidrocarboneto oxigênio     dif   Y.pred     residuo
## 1            0.99    90.01 -2.1505 89.08132  0.92868082
## 2            1.02    89.05 -3.1105 89.52974 -0.47974357
## 3            1.15    91.43 -0.7305 91.47292 -0.04291593
## 4            1.29    93.74  1.5795 93.56556  0.17443691
## 5            1.46    96.73  4.5695 96.10663  0.62336535
## 6            1.36    94.45  2.2895 94.61189 -0.16188668
## 7            0.87    87.59 -4.5705 87.28762  0.30237839
## 8            1.23    91.77 -0.3905 92.66871 -0.89871431
## 9            1.55    99.42  7.2595 97.45191  1.96809217
## 10           1.40    93.65  1.4895 95.20979 -1.55978587
## 11           1.19    93.54  1.3795 92.07082  1.46918488
## 12           1.15    92.52  0.3595 91.47292  1.04708407
## 13           0.98    90.56 -1.6005 88.93184  1.62815562
## 14           1.01    89.54 -2.6205 89.38027  0.15973123
## 15           1.11    89.85 -2.3105 90.87502 -1.02501674
## 16           1.20    90.39 -1.7705 92.22029 -1.83028992
## 17           1.26    93.25  1.0895 93.11714  0.13286130
## 18           1.32    93.41  1.2495 94.01399 -0.60398749
## 19           1.43    94.98  2.8195 95.65821 -0.67821026
## 20           0.95    87.33 -4.8305 88.48342 -1.15341999

• Intervalos de confiança para os parâmetros do modelo de regressão linear simples

confint(fit)

##                   2.5 %   97.5 %
## (Intercept)    70.93555 77.63108
## hidrocarboneto 12.18107 17.71389

• Intervao de confiança para uma resposta média ao valor de \(x_0\)

new.pureza <- data.frame(hidrocarboneto=1.30)
predict(fit, newdata = new.pureza, interval = 'confidence')

##        fit      lwr      upr
## 1 93.71504 93.12911 94.30097

• Intervao de confiança para uma resposta média a uma futura observação \(x_f\)

new.pureza1 <- data.frame(hidrocarboneto=1.60)
predict(fit, newdata = new.pureza1, interval = 'prediction')

##        fit      lwr      upr
## 1 98.19928 95.60691 100.7917

Comparaçõers de Retas de Regressão

Resolução do Exemplo 5.2.1