Para conseguirmos obter a inversa de uma matriz, podemos usar diferentes comandos do R solve()
, ginv()
, etc, porém, ao tentar utilizar esses comando, ocorre um erro onde achamos que o problema deve ser porque há valores negativos na matriz, já que na descrição dos comandos eles especificam que devemos achar a inversa para matrizes positivo-definidas. Como conseguimos resolver essa questão?
-
Você poderia exemplificar com código, ou com o erro que está encontrando.– Lucas Lima23/03/2014 às 16:45
-
X_i<-chol2inv(X), onde X é uma matriz de dimensão 1000x3 Usamos o comando chol2inv() da nossa matriz original na tentativa de achar o valor de (X'X)^-1 e encontramos uma matriz de dimensão 3x3, que é o esperado. Certo? Dado que esse problema foi resolvido, precisamos criar um vetor y com dimensão 1000x1, y sendo uma função (processo gerador de dados). Podemos definir y como a pgd (y<-a+bx+cz) e em seguida definir suas dimensões com (dim(y)=c(1000,1,1)? Fizemos isso e não temos certeza se está correto.– felipe_Milla23/03/2014 às 18:15
-
Edite sua pergunta, e acrescente lá.– Lucas Lima23/03/2014 às 18:17
-
Felipe, não está muito claro qual é o problema, se você conseguir colocar um exemplo do código fica mais fácil de ajudar. Veja aqui como melhorar a pergunta: meta.pt.stackoverflow.com/questions/824/…– Carlos Cinelli23/03/2014 às 18:29
2 Respostas
Felipe, não está muito claro qual é o problema, o ideal seria você colocar o código com o (possível) erro. Mas seguem algumas maneiras de se inverter uma matriz no R
:
Gerando uma base de dados de exemplo:
### constante
c <- rep(1, 1000)
### x ~ N(1,2)
set.seed(1)
x <-rnorm(1000, 1,2)
### w ~ N(1,2)
set.seed(2)
w <- rnorm(1000, 1, 2)
X = cbind(c,x,w)
Invertendo (X'X):
### com solve:
a <- solve(t(X)%*%X)
### com ginv:
b <-ginv(t(X)%*%X)
### com choleski:
c <- chol2inv(chol(t(X)%*%X))
### com qr:
d <-qr.solve(qr(t(X)%*%X))
Comparando todas para ver que são iguais:
### primeiro vamos tirar os dimnames para comparar somente os numeros
dimnames(a) <- NULL
dimnames(d) <- NULL
### comparando
all.equal(a,b,c,d)
[1] TRUE
É preciso usar o comando all.equal
por causa de erros de ponto flutuante.
No caso de estimar os coeficientes de uma regressão (se você não quer usar a função lm
), você pode usar qr.coef
:
### y = 10 + 5x + 2w + u, u~N(0,1)
set.seed(3)
y <- 10 + 5*x +2*w + rnorm(1000)
### estimando os betas
B <- qr.coef(qr(X), y)
B
c x w
9.994012 5.019442 1.994124
se você que inverter X'X
, a maneira mais simples é usando solve()
:
# crio uma matriz qualquer (3 colunas, 100 linhas)
X = cbind(1, rnorm(100), runif(100))
# invertenddo X'X
a = solve(t(X)%*%X)
# invertendo a decomposição QR
b = chol2inv(X)
# é a mesma coisa?
a==b
[,1] [,2] [,3]
[1,] FALSE FALSE FALSE
[2,] FALSE FALSE FALSE
[3,] FALSE FALSE FALSE