Estatística Computacional II - Lista 2

Aluno: Rafael Cabral Fernandez

Introdução

Este documento refere-se a lista 2 de exercícios proposta pelo professor Gustavo Rocha, da cadeira de Estatística Computacional 2, como cômputo de ensino da disciplina.

a)

Considerando o método Gauss-Seidel, utilize 15 iterações e registre os valores obtidos para as soluções das iterações do sistema:

\(\begin{bmatrix} 2 & 1\\ 5 & 7 \end{bmatrix}x\) = \(\begin{bmatrix} 11\\ 13 \end{bmatrix}\)

Considere o vetor nulo como ponto inicial.

Início:

A <- matrix(c(2,5,1,7),2,2)
#Matriz de coeficientes

b <- matrix(c(11,13),2,1)
#Vetor de resposta

xinicial <- matrix(c(0,0),2,1)
#Matriz de solução inicial

L <- matrix(c(0,5,0,0),2,2)
D <- rbind(c(2,0),c(0,7))
U <- matrix(c(0,0,1,0),2,2)
#Matriz decomposta

B <- -1*solve(L+D)
C <- solve(L+D)
#Componentes da Iteração


vetsol <- list()
#Uma lista que contenha os vetores soluções


vetsol[[1]] <- B%*%U%*%xinicial + C%*%b
#Vamos colocar aqui a primeira iteração


for(i in 2:16){
  
  vetsol[[i]] <- B%*%U%*%vetsol[[i-1]] + C%*%b
  print(vetsol[[i]])
  #Equação Interativa
}

##           [,1]
## [1,]  6.535714
## [2,] -2.811224
##           [,1]
## [1,]  6.905612
## [2,] -3.075437
##           [,1]
## [1,]  7.037719
## [2,] -3.169799
##         [,1]
## [1,]  7.0849
## [2,] -3.2035
##           [,1]
## [1,]  7.101750
## [2,] -3.215536
##           [,1]
## [1,]  7.107768
## [2,] -3.219834
##           [,1]
## [1,]  7.109917
## [2,] -3.221369
##           [,1]
## [1,]  7.110685
## [2,] -3.221918
##           [,1]
## [1,]  7.110959
## [2,] -3.222113
##           [,1]
## [1,]  7.111057
## [2,] -3.222183
##           [,1]
## [1,]  7.111092
## [2,] -3.222208
##           [,1]
## [1,]  7.111104
## [2,] -3.222217
##           [,1]
## [1,]  7.111109
## [2,] -3.222220
##           [,1]
## [1,]  7.111110
## [2,] -3.222222
##           [,1]
## [1,]  7.111111
## [2,] -3.222222

b)

Seja a matriz:

\[\begin{bmatrix} 1 & 2 & 1\\ 2 & 8 & 4\\ 1 & 4 & 11 \end{bmatrix}\]

1. Encontre seu fator de Cholesky𝑼 e compare com o resultado obtido pela função ‘chol’ do R.

A <- rbind(c(1,2,1),c(2,8,4),c(1,4,11))
#Matriz que desejamos fatorar

U <- matrix(rep(0,9),3,3)
#Por necessidade, devemos criar uma matriz U "vazia"


choleskao <- function(A,U){
  #install.packages("corpcor")
  #library(corpcor)
  #Poderiamos verificar se a matriz é positiva definida através desse pacote
  
  timesup <- proc.time()  
  
  if(isSymmetric.matrix(A)==0) stop("A matriz precisa ser simétrica")
  #if(is.positive.definite(A)==0) stop("A matriz precisa ser positiva definida")

  
  #Iterador, transcrição direta do algoritmo
  for (i in 1:3){
    U[i,i]=sqrt(A[i,i]-sum(U[1:(i-1),i]^2))
    if (i<3){
      for (j in (i+1):3){
        U[i,j]=(A[i,j]-sum(U[1:(i-1),i]*U[1:(i-1),j]))/U[i,i]
        
        timesend <- proc.time() - timesup
        
      }
      
    }
      print(U)
  }

  print(timesend)
}

choleskao(A,U)

##      [,1] [,2] [,3]
## [1,]    1    2    1
## [2,]    0    0    0
## [3,]    0    0    0
##      [,1] [,2] [,3]
## [1,]    1    2    1
## [2,]    0    2    1
## [3,]    0    0    0
##      [,1] [,2] [,3]
## [1,]    1    2    1
## [2,]    0    2    1
## [3,]    0    0    3
##    user  system elapsed 
##       0       0       0

2. Encontre a solução do sistema \(Ax = (11,13)'\) considerando o fator de Cholesky U obtido em a)

A <- rbind(c(1,2,1),c(2,8,4),c(1,4,11))
#Matriz original


U <- chol(A)
#Como dá o mesmo resultado do algoritmo feito em A, usarei a função chol


b <- matrix(c(-1,0,1),3,1)
#Vetor de respostas 


Ulinha <- t(U)
#Transposta do fator de Cholesky


z <- solve(Ulinha,b)
#Resolvendo o primero sistema para Z

soluction <- solve(U,z)
#Resolvendo o segundo sistema para X


soluction

##            [,1]
## [1,] -2.0000000
## [2,]  0.4444444
## [3,]  0.1111111

c)

Considerando o método de Newton-Raphson, encontre as estimativas de máxima verossimilhança dos parâmetros de um modelo linear simples tal que:

\(Y_i = \beta_o + \beta_1x_i + \epsilon_i\), \(\epsilon_i\sim n(0,\sigma^2)\)

Supondo erros normalmente distribuídos e não-correlacionados, ganhando a propriedade da independência e portanto:

\(Y_i\sim n(\beta_o + \beta_1x_i,\sigma^2)\)

Teremos que a função de verossimilhança é dada por:

\(L(\beta_o,\beta_1,\sigma^2|y_i)=(2\pi\sigma^2)^{-n/2}exp{\cfrac{-1}{2\sigma^2}\sum_{i=1}^{n}[y_i-(\beta_o+\beta_1x_i)]^2}\)

Ademais, devemos resolver o sistema que consite em:

\[\begin{cases} \cfrac{\partial L}{\partial \beta_0} = 0 = u_1\\ \cfrac{\partial L}{\partial \beta_1} = 0 = u_2\\ \cfrac{\partial L}{\partial \sigma^2} = 0 = u_3 \end{cases}\]

Obteremos os seguintes resultados:

\(u_1 = \cfrac{266.43-100\beta_o - 160\beta_1}{\sigma^2}\)

\(u_1 = \cfrac{704.34 - 160\beta_o - 354\beta_1}{\sigma^2}\)

\(u_3 = \cfrac{-100}{\sigma^2}+\cfrac{2388.56-532.86\beta_0-1408.78\beta_1+320\beta_0\beta_1+100(\beta_0)^2+354(\beta_1)^2)}{(\sigma^2)^2}\)

Agora, precisamos achar os valores da matriz jacobiana, determinada por:

\[\begin{bmatrix} \cfrac{\partial u_1}{\partial \beta_0} & \cfrac{\partial u_1}{\partial \beta_1} &\cfrac{\partial u_1}{\partial \sigma^2} \\ \cfrac{\partial u_2}{\partial \beta_0} & \cfrac{\partial u_2}{\partial \beta_1} &\cfrac{\partial u_2}{\partial \sigma^2} \\ \cfrac{\partial u_3}{\partial \beta_0} & \cfrac{\partial u_3}{\partial \beta_1} &\cfrac{\partial u_3}{\partial \sigma^2} \end{bmatrix}\]

Sendo, portanto:

\(\cfrac{\partial u_1}{\partial \beta_0} = \cfrac{-100}{\sigma^2}\)

\(\cfrac{\partial u_1}{\partial \beta_1} = \cfrac{-160}{\sigma^2}\)

\(\cfrac{\partial u_1}{\partial \sigma^2} = \cfrac{100\beta_0 - 160\beta_1 - 266.43}{(\sigma^2)^2}\)

\(\cfrac{\partial u_2}{\partial \beta_0} = \cfrac{-100}{\sigma^2}\)

\(\cfrac{\partial u_2}{\partial \beta_1} = \cfrac{-160}{\sigma^2}\)

\(\cfrac{\partial u_2}{\partial \sigma^2} = \cfrac{100\beta_0 - 160\beta_1 - 266.43}{(\sigma^2)^2}\)

\(\cfrac{\partial u_3}{\partial \beta_0} = \cfrac{-532.86+320\beta_1+200\beta_0}{(\sigma^2)^2}\)

\(\cfrac{\partial u_3}{\partial \beta_1} = \cfrac{-1408.78+320\beta_0+708\beta_1}{(\sigma^2)^2}\)

\(\cfrac{\partial u_3}{\partial \sigma^2} = \cfrac{100}{(\sigma^2)^2}-\cfrac{2(2388.56-532.86\beta_0-1408.78\beta_1+320\beta_0\beta_1+100(\beta_0)^2+354(\beta_1)^2)}{(\sigma^2)^3}\)

De posse de todas as funções necessárias, podemos agora formular o programa:

inicial <- matrix(c(0,0,5),3,1)
#vetor de pontos iniciais



#Função de avalia os pontos nas funções u1, u2 e u3
u_eval <- function(theta){
  
u1 <-(266.93-100*theta[1]-160*theta[2])/(theta[3])

u2 <-(704.39-160*theta[1]-354*theta[2])/(theta[3])

u3 <-((-100/(2*theta[3]))+(2388.56-2*theta[1]*266.93-2*theta[2]*704.39+2*theta[1]*theta[2]*160+100*(theta[1]^2)+(theta[2]^2)*354)/(2*(theta[3]^2)))

u  <- matrix(c(u1,u2,u3),3,1) 
  
}


#Função que avalia os pontos no jacobiano
jacob_eval <- function(theta){
  
  jacob <- matrix(rep(0,9),3,3)
  
  linha1 <- c(-100/theta[3],-160/theta[3],(-266.93+100*theta[1]+160*theta[2])/(theta[3]^2))
  linha2 <- c(-160/theta[3],-354/theta[3],(-704.39+160*theta[1]+354*theta[2])/(theta[3]^2))
  linha3 <- c((-266.93+100*theta[1]+160*theta[2])/(theta[3]^2),(-704.39+160*theta[1]+354*theta[2])/(theta[3]^2),(50/(theta[3]^2))-(2388.56-theta[1]*2*266.93-2*theta[2]*704.39+2*theta[1]*theta[2]*160+100*(theta[1]^2)+(theta[2]^2)*354)/(theta[3]^3))
  
  jacob <- rbind(linha1,linha2,linha3)

}



library(MASS)
#Se não usar, da pau


vetsol <- list()
#Uma lista que contenha os vetores soluções


vetsol[[1]] <- inicial - ginv(jacob_eval(inicial))%*%u_eval(inicial)
#Vamos colocar aqui a primeira iteração



for(i in 2:16){
  vetsol[[i]] <- vetsol[[i-1]] - ginv(jacob_eval(vetsol[[i-1]]))%*%u_eval(vetsol[[i-1]])
  print(vetsol[[i]])
  #Equação Interativa
}

##            [,1]
## [1,] -1.5869148
## [2,]  2.4166637
## [3,]  0.5730682
##            [,1]
## [1,] -1.7320113
## [2,]  2.6376268
## [3,]  0.8394955
##           [,1]
## [1,] -1.798672
## [2,]  2.739142
## [3,]  1.235095
##           [,1]
## [1,] -1.830642
## [2,]  2.787828
## [3,]  1.805535
##           [,1]
## [1,] -1.845912
## [2,]  2.811083
## [3,]  2.606203
##           [,1]
## [1,] -1.853038
## [2,]  2.821935
## [3,]  3.686079
##           [,1]
## [1,] -1.856216
## [2,]  2.826774
## [3,]  5.048694
##           [,1]
## [1,] -1.857516
## [2,]  2.828754
## [3,]  6.575726
##           [,1]
## [1,] -1.857962
## [2,]  2.829434
## [3,]  7.942145
##           [,1]
## [1,] -1.858068
## [2,]  2.829595
## [3,]  8.722832
##           [,1]
## [1,] -1.858079
## [2,]  2.829612
## [3,]  8.905850
##           [,1]
## [1,] -1.858080
## [2,]  2.829612
## [3,]  8.913852
##           [,1]
## [1,] -1.858080
## [2,]  2.829612
## [3,]  8.913866
##           [,1]
## [1,] -1.858080
## [2,]  2.829612
## [3,]  8.913866
##           [,1]
## [1,] -1.858080
## [2,]  2.829612
## [3,]  8.913866

d)

Encontre a raiz da função \(f(x) = 2x^2 - 5\) para \(x>0\), considerando o método de Newton-Raphson e o ponto inicial \(x^{(o)} = 1\). Execute o algoritmo enquanto \(|x^{j}-x^{k-1}| > 0,0001\). Registre os valores obtidos para \(x^{i}\) e as diferenças absolutas \(|x^{j}-x^{k-1}|\)

u <- function(x){
  2*(x^2) - 5
  }
#Função u = f(x)


ulinha <- function(x){
  4*x
}
#Derivada de u = f(x)


x0 <- 1
#Ponto inicial


erro <- vector()
#Vetor para alocar os erros


xj <- vector()
#Vetor para alocar as raízes


xj[1]<- x0 - u(x0)/ulinha(x0)
xj[2]<- xj[1] - u(xj[1])/ulinha(xj[1])
#Os 2 primeiros valores são imputadores manualmente para que o while funcione
#Caso contrário, ele não terá valores para avaliar


i <- 2

#Condicional imposta pela questão
while (abs(xj[i]-xj[(i-1)])>0.0001) {
  xj[i+1] <- xj[i] - u(xj[i])/ulinha(xj[i]) 
  i <- i + 1
}


#Iterador para o calculo dos erros
for (j in 2:(length(xj)+1)){
  erro[j] <- abs(xj[j]-xj[j-1])  
  j <- j+1
}


#Criando um data frame para alocar erros e respostas
tab <-as.data.frame(cbind(xj,erro))
tab[is.na(tab)] <- 0
tab = tab[-5,]


#Usando o Kable para tornar a resposta mais apresentável
library(knitr)
kable(tab[1:dim(tab)[1], ],)

xj	erro
1.750000	0.0000000
1.589286	0.1607143
1.581160	0.0081260
1.581139	0.0000209