Aula de Modelos Lineares - Unidade II - continuação
1 Pacotes Utilizados na disciplina
require(matrixcalc)
require(Matrix)
require(ibd)
require(pracma)
require(matlib)
require(MASS)
require(tidyverse)2 Solução de mínimos quadrados
Sistemas lineares são comumente usados como modelos matemáticos em diversas áreas. No entanto, experimentos aleatórios ou erros de medição resultam em um sistema inconsistente, ou seja, sem soluções exatas. Nessas situações, o Método dos Mínimos Quadrados oferece uma solução aproximada.
2.1 Solução aproximada (de mínimos quadrados):
A solução de mínimos quadrados é dada por: \[ \begin{array}{l} \hat{\boldsymbol{x}}=\boldsymbol{A}^+.\boldsymbol{b}, \mbox{ou equivalentemente },\\ \hat{\boldsymbol{x}}=(\boldsymbol{A}^t\boldsymbol{A})^{-1}\boldsymbol{A}^t \boldsymbol{b}\\ \Rightarrow \mbox{ as soluções são idênticas sempre que } (\boldsymbol{A}^t \boldsymbol{A})^{-1}\mbox{ existir}\\ \Rightarrow \mbox{ é a solução que minimiza a norma euclidiana}: ||\boldsymbol{A}\boldsymbol{x}-\boldsymbol{b}|| \\ \mbox{ significa que estamos minimizando os quadrados das distâncias (ou erros)}\\ \end{array} \]
Exemplo 1: (Reinprecht, 2022) Seja o sistema de equações abaixo (inconsistente): \[ \begin{array}{l} \left\{ \begin{array}{l} 2 x_1+3 x_2=1\\ 3x_1+5x_2=-2\\ -2x_1-4x_2=-1 \end{array} \right. \Rightarrow \underbrace{ \begin{pmatrix} 2 & 3\\ 3 & 5\\ -2 & -4 \end{pmatrix}}_{\boldsymbol{A}} \underbrace{ \begin{pmatrix} x_1\\ x_2 \end{pmatrix}}_{\boldsymbol{x}} \underbrace{ \begin{pmatrix} 1\\ -2\\ -1 \end{pmatrix}}_{\boldsymbol{b}}\\ \square \mbox{ Inversa generalizada: } \boldsymbol{A}^{+}=\begin{pmatrix} 13/9 & 5/9 & 16/9\\ -7/9 & -2/9 & -10/9 \end{pmatrix}\\ \square \mbox{ consistência (não há):} \boldsymbol{A}\boldsymbol{A}^{+}\boldsymbol{b}\ne \boldsymbol{b}\\ \square \mbox{ Solução aproximada: } \hat{\boldsymbol{x}}= \begin{pmatrix} -13/9\\ 7/9 \end{pmatrix} \end{array} \]
\[ \begin{array}{l} \square \mbox{ vetor de resíduos para esta solução:}\\ \boldsymbol{\hat{\epsilon}}=\boldsymbol{b}-\boldsymbol{A}.\hat{\boldsymbol{x}} = \begin{pmatrix} 1\\ -2\\ -1 \end{pmatrix}-\begin{pmatrix} 2 & 3\\ 3 & 5\\ -2 & -4 \end{pmatrix} \begin{pmatrix} -13/9 \\ 7/9 \end{pmatrix}\\ =\begin{pmatrix} 1-(-26/9+21/9)\\ -2-(-39/9+35/9)\\ -1-(26/9-28/9) \end{pmatrix}= \begin{pmatrix} 14/9\\ -14/9\\ -7/9 \end{pmatrix}\\ \square \mbox{ soma de quadrados dos resíduos correspondente}:\\ \Rightarrow \sum \limits_{i=1}^3 \hat{\epsilon_i}=(14/9)^2+(14/9)^2+(7/9)^2=\frac{441}{81}\approx 5.4444\\ \mbox{de forma que } ||\boldsymbol{b}-\boldsymbol{A}\boldsymbol{\hat{x}}||= \sqrt{\frac{441}{81}}\approx 2,33 \\ \mbox{ também denominado desvio padrão estimado (no contexto da estatística)} \end{array} \]
- Cálculos auxiliares:
- Sistema linear em forma matricial:
A=matrix(c(2,3,3,5,-2,-4),ncol=2,byrow=T)
b=matrix(c(1,-2,-1),ncol=1)
showEqn(A, b)## 2*x1 + 3*x2 = 1
## 3*x1 + 5*x2 = -2
## -2*x1 - 4*x2 = -1- Não existe inversa:
try(solve(A))## Error in solve.default(A) : 'a' (3 x 2) deve ser quadrada- Obtenção da inversa generalizada de Moore Penrose:
A_mais=fractions(ginv(A))
A_mais## [,1] [,2] [,3]
## [1,] 13/9 5/9 16/9
## [2,] -7/9 -2/9 -10/9fractions(solve(t(A)%*%A)%*%t(A))## [,1] [,2] [,3]
## [1,] 13/9 5/9 16/9
## [2,] -7/9 -2/9 -10/9- Verificação de consistência:
- não há
all.equal(A%*%A_mais%*%b,b)## [1] "Mean relative difference: 3.181818"- Solução aproximada do sistema: \(\hat{\boldsymbol{x}}=\boldsymbol{A}^+ \boldsymbol{b}\)
x_hat=fractions(A_mais%*%b)
x_hat## [,1]
## [1,] -13/9
## [2,] 7/9- que também pode ser obtida pelo solve de sistemas lineares:
fractions(solve(qr(A, LAPACK=TRUE), b))## [,1]
## [1,] -13/9
## [2,] 7/9- Os resíduos são estimados:
e_hat=b-A%*%x_hat
fractions(e_hat)## [,1]
## [1,] 14/9
## [2,] -14/9
## [3,] -7/9- E soma dos quadrados dos resíduos também (ou norma euclidiana)
norm(b-A%*%x_hat,type="F")## [1] 2.333333- Equivalentemente, utilizamos a ótica de modelo linear sem intercepto, obtendo os mesmos resultados:
m1=lm(b ~ 0 + A)
m1##
## Call:
## lm(formula = b ~ 0 + A)
##
## Coefficients:
## A1 A2
## -1.4444 0.7778- E os resíduos e desvio padrão estimado vem dos resultados do modelo linear sem intercepto:
fractions(residuals(m1)) #para visualizar o vetor de resíduos## 1 2 3
## 14/9 -14/9 -7/9anova(m1) #para visualizar a soma de quadrados dos resíduos## Analysis of Variance Table
##
## Response: b
## Df Sum Sq Mean Sq F value Pr(>F)
## A 2 0.5556 0.2778 0.051 0.9526
## Residuals 1 5.4444 5.4444summary(m1) #para visualizar o desvio padrão estimado##
## Call:
## lm(formula = b ~ 0 + A)
##
## Residuals:
## 1 2 3
## 1.5556 -1.5556 -0.7778
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## A1 -1.4444 5.4997 -0.263 0.836
## A2 0.7778 3.2069 0.243 0.849
##
## Residual standard error: 2.333 on 1 degrees of freedom
## Multiple R-squared: 0.09259, Adjusted R-squared: -1.722
## F-statistic: 0.05102 on 2 and 1 DF, p-value: 0.9526Exemplo 2: (Reinprecht, 2022) Aplicação com banco de dados: seja o banco de dados abaixo, ajuste uma reta de regressão de mínimos quadrados com intercepto do tipo: \(f(x)=\beta_1+\beta_2.x\). Vamos atentar para a mudança de nomenclatura na forma matricial do sistema!
dados=data.frame(X=c(1,2,3,4,5),
Y=c(3.8,7.3,10.2,12.7,15.5))
dados## X Y
## 1 1 3.8
## 2 2 7.3
## 3 3 10.2
## 4 4 12.7
## 5 5 15.5Resolucão gráfica:
modelo <- lm(Y ~ X, data = dados)
ggplot(dados, aes(x = X, y = Y)) +
geom_point(aes(color = "Pontos"), alpha = 0.6) + # Gráfico de dispersão
geom_smooth(method = "lm", se = FALSE, aes(color = "Reta Estimada")) + # Reta de regressão
theme(legend.title = element_blank()) + # Remove o título da legenda
geom_text(x = 1.5, y = 15, label = paste("f(x) = ", round(coef(modelo)[1], 2), " + ", round(coef(modelo)[2], 2), "x"), hjust = 0) # Adicionar a fórmula da reta
Resolução detalhada:
- O sistema linear correspondente é: \[ \begin{array}{l} \left\{ \begin{array}{l} \beta_1+\beta_2=3.8\\ \beta_1+2\beta_2=7.3\\ \beta_1+3\beta_2=10.2\\ \beta_1+4\beta_2=12.7\\ \beta_1+5\beta_2=15.5 \end{array} \right. \Rightarrow \underbrace{ \begin{pmatrix} 1 & 1\\ 1 & 2\\ 1 & 3\\ 1 & 4\\ 1 & 5\\ \end{pmatrix}}_{\boldsymbol{X}} \underbrace{ \begin{pmatrix} \beta_1\\ \beta_2 \end{pmatrix}}_{\boldsymbol{\beta}}= \underbrace{ \begin{pmatrix} 3.8\\ 7.3\\ 10.2\\ 12.7\\ 15.5 \end{pmatrix}}_{\boldsymbol{y}}\\ \square \mbox{ Inversa generalizada: } \boldsymbol{A}^{+}=\begin{pmatrix} 4/5 & 1/2 & 1/5 & -1/10 & -2/5\\ -1/5 & -1/10 & 0 & 1/10 & 1/5 \end{pmatrix}\\ \square \mbox{ Solução aproximada: } \hat{\boldsymbol{\beta}}= \begin{pmatrix} 1.26\\ 2.88 \end{pmatrix}\\ \Rightarrow f(x)=1.26+2.88.x \end{array} \]
\[ \begin{array}{l} \square \mbox{ vetor de resíduos: } \boldsymbol{\hat{\epsilon}}= \begin{pmatrix} -17/50 \\ 7/25 \\ 3/10 \\-2/25 \\ -4/25 \end{pmatrix}\\ \square \mbox{ soma de quadrado de resíduos}: \sum \limits_{i=1}^5 \hat{\boldsymbol{\epsilon_i}}=0,316\\ \square \mbox{ desvio padrão estimado : }\hat{\sigma}=0,3246 \end{array} \]
- Cálculos auxiliares:
- Sistema linear em forma matricial:
X=cbind(rep(1,5),dados$X)
Y=matrix(dados$Y,ncol=1)
showEqn(X, Y,c("beta1","beta2"))## 1*beta1 + 1*beta2 = 3.8
## 1*beta1 + 2*beta2 = 7.3
## 1*beta1 + 3*beta2 = 10.2
## 1*beta1 + 4*beta2 = 12.7
## 1*beta1 + 5*beta2 = 15.5- Não existe inversa:
try(solve(X))## Error in solve.default(X) : 'a' (5 x 2) deve ser quadrada- Obtenção da inversa generalizada de Moore Penrose:
X_mais=fractions(ginv(X))
X_mais## [,1] [,2] [,3] [,4] [,5]
## [1,] 4/5 1/2 1/5 -1/10 -2/5
## [2,] -1/5 -1/10 0 1/10 1/5- Solução aproximada do sistema: \(\hat{\boldsymbol{x}}=\boldsymbol{A}^+ \boldsymbol{b}\)
X_mais%*%Y## [,1]
## [1,] 1.26
## [2,] 2.88- ou equivalentemente
fractions(solve(qr(X, LAPACK=TRUE), Y))## [,1]
## [1,] 63/50
## [2,] 72/25- Por fim, a resolução utilizando o comando lm (linear
model):
- pega os valores direto do banco de dados (o banco de dados não tem coluna de un’s)
m2=lm(Y ~ X,data=dados)
m2##
## Call:
## lm(formula = Y ~ X, data = dados)
##
## Coefficients:
## (Intercept) X
## 1.26 2.88- E os resíduos e desvio padrão estimado vem dos resultados do modelo linear com intercepto:
fractions(residuals(m2)) #para visualizar o vetor de resíduos## 1 2 3 4 5
## -17/50 7/25 3/10 -2/25 -4/25anova(m2) #para visualizar a soma de quadrados dos resíduos## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X 1 82.944 82.944 787.44 9.935e-05 ***
## Residuals 3 0.316 0.105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m2) #para visualizar o desvio padrão estimado##
## Call:
## lm(formula = Y ~ X, data = dados)
##
## Residuals:
## 1 2 3 4 5
## -0.34 0.28 0.30 -0.08 -0.16
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.2600 0.3404 3.702 0.0342 *
## X 2.8800 0.1026 28.061 9.93e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3246 on 3 degrees of freedom
## Multiple R-squared: 0.9962, Adjusted R-squared: 0.9949
## F-statistic: 787.4 on 1 and 3 DF, p-value: 9.935e-053 Referências
[1] DE LUNA, J. G., & DE OLINDA, R. A. (2015). Introdução a Modelos Lineares. Livraria da Física, 1ª edição, 164 p.
[2] FIELLER, N. (2016). Basics of Matrix Algebra for Statistics with R. (Chapman & Hall/CRC The R Series), 1ª edição, 244 p..
[3] GRAYBILL, F.A. Introduction to matrices with applications. Belmont, Califórnia, Wadsworth Pub. Cia, 1969.
[4] RAO, C.Radhakrishna ; MITRA, Sujit Kumar. Generalized Inverse of Matrix and its Applicátion.. New York: Edição Eletrônica, Wiley & Sons Inc, 1971.
[5] REINPRECHT, J.. (2022). A inversa generalizada de uma matriz e aplicações. Dissertação de mestrado da Universidade Estadual Júlio de Mesquita Filho, Rio Claro, SP.