1. No category

Cálculos en Regresión
Felipe Osorio
http://www.ies.ucv.cl/fosorio
Instituto de Estadı́stica
Pontificia Universidad Católica de Valparaı́so
Abril 21, 2015
1 / 19
Estimación en regresión lineal
I El método por defecto para llevar a cabo la estimación de parámetros en
regresión lineal es mediante la descomposición QR.
I También podemos usar los siguientes procedimientos:
I
Descomposición Cholesky.
I
Descomposición SVD.
I
Operador Sweep.
I
Método de gradientes conjugados.
2 / 19
Estimación en regresión lineal: Método Cholesky
Podemos usar el método Cholesky (¿por qué?) para obtener estimaciones de los
coeficientes de regresión, mediante considerar el sistema de ecuaciones:
b = GG> β
b = X>Y ,
X>Xβ
y resolvemos los sistemas triangulares:
Gz = X > Y ,
b = z.
G> β
Observaciones:
I El procedimiento requiere “formar” las matrices X > X y X > Y .
I Estos resultados no pueden ser “aprovechados” para obtener s2 , los valores
predichos y/o residuos.
3 / 19
Estimación en regresión lineal: Método Cholesky
Podemos usar el método Cholesky (¿por qué?) para obtener estimaciones de los
coeficientes de regresión, mediante considerar el sistema de ecuaciones:
b = GG> β
b = X>Y ,
X>Xβ
y resolvemos los sistemas triangulares:
Gz = X > Y ,
b = z.
G> β
Observaciones:
I El procedimiento requiere “formar” las matrices X > X y X > Y .
I Estos resultados no pueden ser “aprovechados” para obtener s2 , los valores
predichos y/o residuos.
3 / 19
Estimación en regresión lineal: Método SVD
Considere la descomposición SVD de X como:
X = U DV ,
Rn×p ,
>
donde U ∈
tal que U U = I p . D = diag(δ1 , . . . , δp ) y V ∈ Op . De este
modo, podemos re-escribir el modelo lineal como:
Y = Xβ + = U DV β + .
Sea α = V β, esto lleva al modelo en forma canónica:
Y = U Dα + ,
Cov() = σ 2 I.
4 / 19
Estimación en regresión lineal: Método SVD
Sea
Z = U >Y ,
Cov(U > ) = σ 2 I.
Entonces, la estimación de parámetros en el modelo canónico es dada por:
b = D −1 z.
α
Finalmente, para obtener el estimador LS en el modelo original, calculamos:
b = V >α
b = V > D −1 U > Y .
β
Observaciones:
I El procedimiento permite la estimación en problemas de rango deficiente,
b = D − z con D − alguna inversa generalizada de D.
mediante α
I Podemos reescribir de forma simple algunos resultados, como:
b = U U > Y = HY ,
b = Xβ
Y
En particular, hii =
u>
i ui ,
b
e = Y − U D α.
i = 1, . . . , n.
5 / 19
Estimación en regresión lineal: Método SVD
Sea
Z = U >Y ,
Cov(U > ) = σ 2 I.
Entonces, la estimación de parámetros en el modelo canónico es dada por:
b = D −1 z.
α
Finalmente, para obtener el estimador LS en el modelo original, calculamos:
b = V >α
b = V > D −1 U > Y .
β
Observaciones:
I El procedimiento permite la estimación en problemas de rango deficiente,
b = D − z con D − alguna inversa generalizada de D.
mediante α
I Podemos reescribir de forma simple algunos resultados, como:
b = U U > Y = HY ,
b = Xβ
Y
En particular, hii =
u>
i ui ,
b
e = Y − U D α.
i = 1, . . . , n.
5 / 19
Alternativas a mı́nimos cuadrados
I Soluciones regularizadas: Regresión ridge.
I Estimación vı́a IRLS:
I
Modelos lineales generalizados.
I
Estimación L1 .
I
Estimación M .
I
Regresión lineal considerando distribuciones con colas pesadas.
6 / 19
Colinealidad en regresión lineal
Considere el modelo de regresión lineal
Y = Xβ + ,
con X ∈ Rn×p , E() = 0 y Cov() = σ 2 I n .
Es bien conocido que cuando X es mal-condicionada, el sistema de ecuaciones
b = X>Y ,
X>Xβ
puede ser inestable (ver, Stewart, 1987 y Belsley, 1991).
7 / 19
Detección de colinealidad en regresión lineal
Considere la descomposición valor singular (SVD) de X,
X = U DV > ,
donde U ∈ Rn×r , tal que U > U = I r , D = diag(δ1 , . . . , δr ) con δ1 ≥ · · · ≥ δr > 0,
V ∈ Rr×r es matriz ortogonal y r = rg(X).
La detección de colinealidad en el modelo lineal puede ser llevada a cabo por medio
del número condición
δ1
,
κ(X) = kXkkX + k =
δr
y κ(X) “grande” es indicador de colinealidad.
8 / 19
Número condición
Considere la matriz
A=
1.000
0.667
0.500
,
0.333
A−1 =
−666
1344
1000
.
−2000
El número condición se define como κ(A) = kAkkA−1 k para k · k alguna norma
matricial.
Por ejemplo1 ,
κ1 (A) = kAk1 kA−1 k1 = (1.667)(3000) = 5001
κ∞ (A) = kAk∞ kA−1 k∞ = (1.500)(3344) = 5016
λ
maxx6=0 kAxk/kxk
1.333375
max κ2 (A) =
=
= 3555.778
=
λmin minx6=0 kAxk/kxk
0.000375
1 kAk
p
= maxkxkp =1 kAxkp , kAk2 =
p
ρ(A> A)
9 / 19
Cemento Portland (Woods, Steinour y Starke, 1932)
Estudio experimental relacionando la emisión de calor durante la producción y
endurecimiento de 13 muestras de cementos Portland. Woods et al. (1932)
consideraron cuatro compuestos para los clinkers desde los que se produce el cemento.
La respuesta (Y ) es la emisión de calor después de 180 dı́as de curado, medido en
calorı́as por gramo de cemento. Los regresores son los porcentajes de los cuatro
compuestos: aluminato tricálcico (X1 ), silicato tricálcico (X2 ), ferrito aluminato
tetracálcico (X3 ) y silicato dicálcico (X4 ).
10 / 19
Cemento Portland (Woods, Steinour y Starke, 1932)
Siguiendo a Woods et al. (1932) consideramos un modelo lineal sin intercepto
(modelo homogéneo). El número condición escalado es κ(X) = 9.432, esto es X es
bien condicionada. (variables centradas, κ(X̃) = 37.106)
Por otro lado, Hald (1952), Gorman y Toman (1966) y Daniel y Wood (1980) adoptan
un modelo con intercepto (modelo no homogéneo). En cuyo caso, κ(X) = 249.578,
sugiriendo la presencia de colinealidad.
El aumento en el número condición se debe a que existe una relación lineal
aproximada, pues
x1 + x2 + x3 + x4 ≈ 100.
de modo que incluir el intercepto causa una colinealidad severa.
11 / 19
Tratamiento de colinealidad en regresión lineal
El estimador ridge (Hoerl y Kennard, 1970),
b = (X > X + λI)−1 X > Y ,
β
λ
λ > 0.
puede ser visto como la solución del problema regularizado (restringido):
min kY − Xβk2 −
β
λ
kβk2 ,
2
o bien como el problema mı́nimos cuadrados con datos aumentados:
X
Y
= √
β+
.
0
∗
λI p
En este contexto λ es un parámetro de regularización (parámetro ridge).
12 / 19
Cemento Portland
Resultados de estimación:
Estimador
LS (No homog.)
LS (homogéneo)
Ridge
β0
62.405
8.587
β1
2.193
1.551
2.105
β2
1.153
0.510
1.065
β3
0.759
0.102
0.668
β4
0.486
-0.144
0.400
σ2
4.047
3.682
4.000
b = ps2 /kbk2 = 0.00768 (Hoerl, Kennard y Baldwin, 1975).
Se utilizó λ
13 / 19
Regresión L1
Uno de los primeros procedimientos robustos en regresión2 corresponde al problema:
min
β
n
X
|Yi − x>
i β|.
i=1
Mı́nimo desvı́o absoluto (LAD) o regresión L1 puede ser planteado como un problema
de programación lineal, considerando las partes positivas y negativas de los residuos,
e+ y e− , respectivamente, y análogamente para β + , β − .
Ası́, el problema puede ser expresado como (Charnes, Cooper y Ferguson, 1955):
min 1> (e+ + e− ),
β
sujeto a: Y = X(β + − β − ) + (e+ − e− ),
con β + , β − , e+ , e− deben ser todos ≥ 0.
Observación: Barrodale y Roberts (1973, 1974) presentan un algoritmo de propósito
especial para resolver éste problema modificando el método simplex y la estructura de
datos requerida.
2 Este método es, de hecho, anterior a LS!
14 / 19
Regresión L1
Uno de los primeros procedimientos robustos en regresión2 corresponde al problema:
min
β
n
X
|Yi − x>
i β|.
i=1
Mı́nimo desvı́o absoluto (LAD) o regresión L1 puede ser planteado como un problema
de programación lineal, considerando las partes positivas y negativas de los residuos,
e+ y e− , respectivamente, y análogamente para β + , β − .
Ası́, el problema puede ser expresado como (Charnes, Cooper y Ferguson, 1955):
min 1> (e+ + e− ),
β
sujeto a: Y = X(β + − β − ) + (e+ − e− ),
con β + , β − , e+ , e− deben ser todos ≥ 0.
Observación: Barrodale y Roberts (1973, 1974) presentan un algoritmo de propósito
especial para resolver éste problema modificando el método simplex y la estructura de
datos requerida.
2 Este método es, de hecho, anterior a LS!
14 / 19
Llamadas telefónicas en Bélgica 1950-73 (Rousseeuw y Leroy, 1987)
200
●
●
150
●
●
●
calls
50
100
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
0
●
●
50
55
60
65
70
year
Ajustes: LS, normal contaminada ( = .15, γ = .10), Student-t (ν = 2.5), L1 .
15 / 19
Regresión L1
I Schlossmacher (1973) originalmente propuso calcular estimadores LAD usando
IRLS.
I Lange y Sinsheimer (1993) y Phillips (2002) identificaron que éste procedimiento
IRLS corresponde a un algoritmo EM.
I Sin embargo, también ha sido reportado que este algoritmo puede ser incapaz de
detectar la observaciones básicas de manera eficiente.
16 / 19
Regresión lineal: función LAD de L1pack
Considere el modelo
√
−1
i ,
Yi = x>
i β + ( 2τi )
i = 1, . . . , n
ind
donde ∼ N (0, σ 2 ) y τi tiene función de densidad
g(τi ) = τi−3 exp{− 21 τi−2 }.
El algoritmo EM procede a llevar a cabo la estimación de β y σ 2 iterativamente
mediante maximizar la función:
Q(θ|θ (k) ) = −
=−
(k)
n
n
1 X (k)
2
log σ 2 −
W (Yi − x>
i β)
2
2φ i=1 i
n
1
log σ 2 −
(Y − Xβ)> W (k) (Y − Xβ)
2
2φ
(k)
donde W (k) = diag(W1 , . . . , Wn ) y los pesos son dados por:
√
(k)
(k)
|,
Wi = E(τi2 |Yi , θ (k) ) = σ (k) / 2|Yi − x>
i β
(k)
para |Yi − x>
| > 0.
i β
Regresión lineal: función LAD de L1pack
Considere el modelo
√
−1
i ,
Yi = x>
i β + ( 2τi )
i = 1, . . . , n
ind
donde ∼ N (0, σ 2 ) y τi tiene función de densidad
g(τi ) = τi−3 exp{− 21 τi−2 }.
El algoritmo EM procede a llevar a cabo la estimación de β y σ 2 iterativamente
mediante maximizar la función:
Q(θ|θ (k) ) = −
=−
(k)
n
1 X (k)
n
2
log σ 2 −
W (Yi − x>
i β)
2
2φ i=1 i
n
1
log σ 2 −
(Y − Xβ)> W (k) (Y − Xβ)
2
2φ
(k)
donde W (k) = diag(W1 , . . . , Wn ) y los pesos son dados por:
√
(k)
(k)
|,
Wi = E(τi2 |Yi , θ (k) ) = σ (k) / 2|Yi − x>
i β
(k)
para |Yi − x>
| > 0.
i β
Algoritmo IRLS
Paso de coeficientes:
1/2
I Calcular r (k) = Y − Xβ (k) y W (k)
(k)
(k)
= diag(W1 , . . . , Wn )
I Obtener δ (k) como solución del problema WLS
1/2
min kW (k)
δ
(r (k) − Xδ)k2
I Actualizar β (k+1) = β (k) + δ (k) .
Paso de escala:
φ(k+1) =
1/2
1
kW (k) r (k+1) k2
n
Criterio de convergencia: basado en el criterio usado en la función glm.fit().
18 / 19
Detalles de la implementación
c 1/2 X, Y ∗ = W
c 1/2 Y y calcular la descomposición QR de X ∗
Sea X ∗ = W
(DGEQRF)
X∗ = Q
R
,
0
Q ∈ On y R ∈ R
p×p
triangular superior,
considere c = Q> Y ∗ , entonces
(DORMQR)
>
> c 1/2
Q e∗ = Q W
(Y − Xβ) =
c1 − Rβ
c2
=
r1
r2
de este modo, δ es solución del sistema triangular
(DTRTRS)
Rδ = r 1
⇒
Rβ
(k+1)
= c1 ,
actualizar β (k+1) = β (k) + δ (DAXPY) y φ(k+1) = kr 2 k2 /n (DNRM2). Finalmente, note que
(DORMQR)
b ∗ = X ∗ β (k+1) = Q
Y
Rβ (k+1)
0
=Q
c1
0