\[ A = QR\] \[Q : \textrm{Matriz ortonormal}\] \[R : \textrm{Matriz triangular superior}\] Ejemplo: \[A= \begin{pmatrix} 12 & -51 & 4\\ 6 & 167 & -68\\ -4 & 24 & -41 \end{pmatrix}\]
Aplicando la descomposición QR, la matriz A se expresa como el producto de dos matrices Q: ortonormal y R: triangular superior
\[ \begin{pmatrix} 12 & -51 & 4\\ 6 & 167 & -68\\ -4 & 24 & -41 \end{pmatrix} = \begin{pmatrix} 6/7 & -69/175& -58/175\\ 3/7 & 158/175 & 6/175\\ -2/7 & 6/35 & -33/35 \end{pmatrix}\begin{pmatrix} 14 & 21 & -14\\ 0 & 175 & -70\\ 0 & 0 & 35 \end{pmatrix}\]
Donde
\[ Q = \begin{pmatrix} 6/7 & -69/175& -58/175\\ 3/7 & 158/175 & 6/175\\ -2/7 & 6/35 & -33/35 \end{pmatrix}\]
\[R = \begin{pmatrix} 14 & 21 & -14\\ 0 & 175 & -70\\ 0 & 0 & 35 \end{pmatrix}\]
Referencia : Calculadora de factorización QR
Fuente: Computing OLS Coefficients Using QR Decomposition
# Funciones para descomposicion QA
inner_prod = function(v1, v2) {
stopifnot(length(v1) == length(v2))
len = length(v1)
res = vector("numeric", length = len)
for (i in 1:len) {
res[i] = v1[i] * v2[i]
}
return(sum(res))
}
inner = function(v1, v2) {
stopifnot(length(v1) == length(v2))
return(sum(v1 * v2))
}
# Thin QR factorization and modified gram-schmidt
qr_gs = function(x) {
n = nrow(x)
p = ncol(x)
q1 = matrix(0, nrow = n, ncol = p)
r1 = matrix(0, nrow = p, ncol = p)
u = matrix(0, nrow = n, ncol = p)
u[, 1] = x[, 1]
for (k in 2:ncol(x)) {
u[, k] = x[, k]
# successive orthogonalization
for (ctr in seq(1, (k - 1), 1)) {
# classical GS
#u[, k] = u[, k] - ((inner(u[, ctr], x[, k]) / inner(u[, ctr], u[, ctr])) * (u[, ctr]))
# modified (stable) GS
u[, k] = u[, k] - ((inner(u[, ctr], u[, k]) / inner(u[, ctr], u[, ctr])) * (u[, ctr]))
}
}
# dividing each column by respective vector (l2) norm (rescaling each vector to have length 1)
q1 = apply(u, 2, function(x) { x / sqrt(inner(x, x)) })
r1 = crossprod(q1, x) # t(q1) %*% x
return(list(q = q1, r = r1, u = u))
}Creamos la matriz A del ejemplo
A = matrix(c(12,-51,4,6,167,-68,-4,24,-41),nrow = 3, ncol = 3, byrow = TRUE)
A## [,1] [,2] [,3]
## [1,] 12 -51 4
## [2,] 6 167 -68
## [3,] -4 24 -41Procedemos a la decomposición QR utilizando las funciones
Q = qr_gs(A)$q
Q## [,1] [,2] [,3]
## [1,] 0.8571429 -0.3942857 -0.33142857
## [2,] 0.4285714 0.9028571 0.03428571
## [3,] -0.2857143 0.1714286 -0.94285714R = qr_gs(A)$r
R## [,1] [,2] [,3]
## [1,] 1.400000e+01 2.100000e+01 -14
## [2,] -6.661338e-16 1.750000e+02 -70
## [3,] 0.000000e+00 1.421085e-14 35Comprobamos si QR = A
#Comprobando
A1 = Q%*%R
A1## [,1] [,2] [,3]
## [1,] 12 -51 4
## [2,] 6 167 -68
## [3,] -4 24 -41Sea la matriz :
\[A= \begin{pmatrix} 12 & -51 & 4\\ 6 & 167 & -68\\ -4 & 24 & -41 \end{pmatrix}\]
A = matrix(c(12,-51,4,6,167,-68,-4,24,-41),nrow = 3, ncol = 3, byrow = TRUE)
A## [,1] [,2] [,3]
## [1,] 12 -51 4
## [2,] 6 167 -68
## [3,] -4 24 -41La matriz Q
QR = qr(A)
Q = qr.Q(QR,complete = FALSE)
Q## [,1] [,2] [,3]
## [1,] -0.8571429 0.3942857 0.33142857
## [2,] -0.4285714 -0.9028571 -0.03428571
## [3,] 0.2857143 -0.1714286 0.94285714La matriz R
R = qr.R(QR,complete = FALSE)
R## [,1] [,2] [,3]
## [1,] -14 -21 14
## [2,] 0 -175 70
## [3,] 0 0 -35Comprobamos si QR = A
#Comprobando
A2 = Q%*%R
A2## [,1] [,2] [,3]
## [1,] 12 -51 4
## [2,] 6 167 -68
## [3,] -4 24 -41Verificando que la función que implementa el método Gram-Schmidt es igual a la función de la librería R
all.equal(abs(qr.Q(qr(A),complete = FALSE)) ,abs(qr_gs(A)$q ))## [1] TRUEFunción generadora de los datos para el modelo:
\[Y \sim N (\mu,\sigma^{2}I)\] \[ \mu = XB\] Referencia : Simbolos LaTex
f_ProcesoGeneradorDatos = function(pN,pP,pB1,pB2,pB3,pB4,pSigma){
X = matrix( ,pN,pP)
Y = matrix( ,pN,1)
B = matrix( c(pB1,pB2,pB3,pB4), pP,1)
Xi= matrix( ,1,pP)
for ( i in 1:pN) {
# Simulando las observaciones
Xi[1,1] = 1
Xi[1,2] = runif(1,20,80)
Xi[1,3] = rexp(1,0.01)
Xi[1,4] = rnorm(1,mean=40,sd=5)
X[i,] = Xi
# Y ~ N(Ui,s^2) con Ui = (Xi)*B
Y[i,] = rnorm(n = 1, mean = X[i,]%*%B, sd= pSigma)
}
M = cbind(Y,X)
colnames(M) = c("Y","X1","X2","X3","X4")
return(M)
} # Parámetros poblacionales (no observables)
B1 = 100
B2 = 3.2
B3 = 6.0
B4 = 1.2
Sigma = 8# n observaciones y p parámetros a estimar
n = 50
p = 4
# Datos Observables
M = f_ProcesoGeneradorDatos(n,p,B1,B2,B3,B4,Sigma)
Y = matrix(M[,1],n,1)
colnames(Y) = c("Y")
X = matrix(M[,2:5],n,p)
colnames(X) = c("X1","X2","X3","X4")
head(M,20)## Y X1 X2 X3 X4
## [1,] 639.8637 1 79.06047 40.514269 37.67970
## [2,] 1500.6843 1 47.95932 200.302322 42.57999
## [3,] 1884.1731 1 74.06981 252.548098 35.69676
## [4,] 1536.9498 1 51.65512 204.285501 27.88576
## [5,] 1281.9427 1 40.13379 166.150986 36.56941
## [6,] 1397.1152 1 21.50122 196.891346 46.13737
## [7,] 1368.1266 1 49.97464 177.925326 36.97609
## [8,] 2049.2079 1 50.27072 286.861323 47.30989
## [9,] 1489.6447 1 58.80136 191.980770 40.60992
## [10,] 501.5841 1 34.21762 40.680068 37.01987
## [11,] 433.4451 1 45.54520 22.684612 38.87516
## [12,] 1002.6969 1 21.83937 131.273833 39.46614
## [13,] 493.6100 1 78.39045 16.068281 36.48095
## [14,] 327.9978 1 54.75728 1.163614 38.01646
## [15,] 616.4412 1 65.90679 42.822490 36.15725
## [16,] 916.4486 1 43.37021 105.223347 37.39406
## [17,] 837.0307 1 21.98681 102.461175 38.94412
## [18,] 452.9892 1 61.79955 15.310521 44.64591
## [19,] 2385.4934 1 52.75665 343.243452 43.53826
## [20,] 519.1281 1 22.85887 49.176017 43.70381\[ X = QR^{*} \] \[ Q = \left[Qf \mid Qg \right]\] \[ R^{*} = \left[ \frac{R}{0} \right]\] \[ Q'Y = \left[ \frac{f}{r} \right]\] \[ \hat{\beta} = R^{-1}f\] \[ SCR = r'r \]
# Estimacion por DESCOMPOSICION QR
# Descomposición QR
QR = qr(X)
Q = qr.Q(QR,complete = TRUE)
Ra = qr.R(QR,complete = TRUE)
Qf = qr.Q(QR,complete = FALSE)
R = qr.R(QR,complete = FALSE)
f = matrix((t(Q)%*%Y)[1:p,1],p,1)
r = matrix((t(Q)%*%Y)[(p+1):n,1],n-p,1)
# --------------------------------------------
B_hat = solve(R)%*%f
SCR = t(r)%*%r
SCT = t(Y)%*%Y- n*(mean(Y)^2)
B1_qr = B_hat[1,1]
B2_qr = B_hat[2,1]
B3_qr = B_hat[3,1]
B4_qr = B_hat[4,1]
s_qr = sqrt(SCR/(n-p))
R2_qr = 1 - SCR/SCTLos resultados de estimación por el método de Descomposición QR:
\(\hat{\beta}_{1}\) = 97.6436711
\(\hat{\beta}_{2}\) = 3.2639977
\(\hat{\beta}_{3}\) = 6.006291
\(\hat{\beta}_{4}\) = 1.1655898
\(\hat{\sigma}\) = 7.1279157
\(R^{2}\) = 0.9998669
\[ \hat{\beta} = (X'X)^{-1}X'Y \] \[ \hat{Y} = X\hat{\beta}\] \[ \hat{e} = Y - \hat{Y}\] \[ SCR = \hat{e}'\hat{e}\]
# Estimation por MINIMOS CUADRADOS ORDINARIOS
B_hat = solve(t(X)%*%X)%*%(t(X)%*%Y)
Y_hat = X%*%B_hat
e_hat = Y - Y_hat
SCR = t(e_hat)%*%e_hat
SCT = t(Y)%*%Y- n*(mean(Y)^2)
B1_ols = B_hat[1,1]
B2_ols = B_hat[2,1]
B3_ols = B_hat[3,1]
B4_ols = B_hat[4,1]
s_ols = sqrt(SCR/(n-p))
R2_ols = 1 - SCR/SCTLos resultados de estimación por el método de Mínimos Cuadrados Ordinarios:
\(\hat{\beta}_{1}\) = 97.6436711
\(\hat{\beta}_{2}\) = 3.2639977
\(\hat{\beta}_{3}\) = 6.006291
\(\hat{\beta}_{4}\) = 1.1655898
\(\hat{\sigma}\) = 7.1279157
\(R^{2}\) = 0.9998669
mrl = lm(Y~X2+X3+X4,data = as.data.frame(M))
B1_lm = coefficients(mrl)[1]
B2_lm = coefficients(mrl)[2]
B3_lm = coefficients(mrl)[3]
B4_lm = coefficients(mrl)[4]
s_lm = sigma(mrl)
R2_lm = summary(mrl)$r.squared\(\hat{\beta}_{1}\) = 97.6436711
\(\hat{\beta}_{2}\) = 3.2639977
\(\hat{\beta}_{3}\) = 6.006291
\(\hat{\beta}_{4}\) = 1.1655898
\(\hat{\sigma}\) = 7.1279157
\(R^{2}\) = 0.9998669
# Cuadro comparativo de resultados
c_Parametro = c('B1','B2','B3','B4','s','R2')
c_Poblacion = c( B1,B2,B3,B4,Sigma, 1.0)
c_qr = c( B1_qr,B2_qr,B3_qr,B4_qr,s_qr,R2_qr)
c_ols = c( B1_ols,B2_ols,B3_qr,B4_ols,s_qr,R2_ols)
c_lm = c( B1_lm,B2_lm,B3_lm,B4_lm,s_qr,R2_lm)
df = data.frame( c_Parametro, c_Poblacion, c_qr, c_ols, c_lm)
title =
"Comparación de los parametros de la Población con los Métodos de Estimación"
knitr::kable(df, col.names = c("Parámetro","Población", "QR","OLS","lm"),
aling = "ccccc",
caption =title )| Parámetro | Población | QR | OLS | lm |
|---|---|---|---|---|
| B1 | 100.0 | 97.6436711 | 97.6436711 | 97.6436711 |
| B2 | 3.2 | 3.2639977 | 3.2639977 | 3.2639977 |
| B3 | 6.0 | 6.0062910 | 6.0062910 | 6.0062910 |
| B4 | 1.2 | 1.1655898 | 1.1655898 | 1.1655898 |
| s | 8.0 | 7.1279157 | 7.1279157 | 7.1279157 |
| R2 | 1.0 | 0.9998669 | 0.9998669 | 0.9998669 |