Considere el siguiente código R:
example <- function(n) {
X <- 1:n
Y <- rep(1,n)
return(lm(Y~X))
}
#(2.13.0, i386-pc-mingw32)
summary(example(7)) #R^2 = .1963
summary(example(62)) #R^2 = .4529
summary(example(4540)) #R^2 = .7832
summary(example(104))) #R^2 = 0
#I did a search for n 6:10000, the result for R^2 is NaN for
#n = 2, 4, 16, 64, 256, 1024, 2085 (not a typo), 4096, 6175 (not a typo), and 8340 (not a typo)
Mirar http://svn.r-project.org/R/trunk/src/appl/dqrls.f ) no me ayudó a entender lo que está sucediendo, porque no conozco Fortran. En otra pregunta , se respondió que los errores de tolerancia de máquina de coma flotante eran los culpables de los coeficientes para X que están cerca, pero no del todo 0.
es mayor cuando el valor deestá más cerca de 0. Pero ...coef(example(n))["X"]
- ¿Por qué hay un valor ?
- ¿Qué (específicamente) lo está determinando?
- ¿Por qué la progresión aparentemente ordenada de los
NaN
resultados? - ¿Por qué las violaciones de esa progresión?
- ¿Qué de esto es el comportamiento 'esperado'?
r
regression
russellpierce
fuente
fuente
Y <- rep(1,n)+runif(n)*ynoise
), eso sería interesante :-)apply(as.matrix(2:17), 1, function(n){example(n)$coefficients[-1]})
. (Mis resultados, en un Win 7 x64 Xeon, varían de -8e-17 a + 3e-16; aproximadamente la mitad son ceros verdaderos). Por cierto, la fuente Fortran no es de ayuda: es solo un contenedor para dqrdc; ese es el código que quieres mirar.Respuestas:
Como dice Ben Bolker, la respuesta a esta pregunta se puede encontrar en el código para
summary.lm()
.Aquí está el encabezado:
Entonces,
x <- 1:1000; y <- rep(1,1000); z <- lm(y ~ x)
eche un vistazo a este extracto ligeramente modificado:Para responder una pregunta con una pregunta: ¿qué extraemos de esto? :)
mss
rss
mss
rss
0/0
NaN
2^(1:k)
Actualización 1: Aquí hay un buen hilo de R-help que aborda algunas de las razones por las cuales las advertencias de subflujo no se abordan en R.
Además, este SO Q&A tiene una serie de publicaciones interesantes y enlaces útiles sobre flujo inferior, aritmética de mayor precisión, etc.
fuente
Tengo curiosidad acerca de su motivación para hacer la pregunta. No puedo pensar en una razón práctica por la que este comportamiento debería importar; La curiosidad intelectual es una razón alternativa (y la OMI mucho más sensata). Creo que no necesita comprender FORTRAN para responder esta pregunta, pero creo que necesita saber sobre la descomposición de QR y su uso en la regresión lineal. Si lo trata
dqrls
como un cuadro negro que calcula una descomposición QR y devuelve información diversa al respecto, entonces puede seguir los pasos ... o simplemente ir directamentesummary.lm
y rastrear para ver cómo se calcula el R ^ 2. En particular:Luego debe volver a
lm.fit
ver que los valores ajustados se calculan comor1 <- y - z$residuals
(es decir, como la respuesta menos los residuos). Ahora puede averiguar qué determina el valor de los residuos y si el valor menos su media es exactamente cero o no, y a partir de ahí averiguar los resultados computacionales ...fuente
mss
rss
fuente