Don’t be loopy! (II)

23 septiembre, 2011 Carlos J. Gil Bellosta 12 comentarios

Continúo en esta la primera de las entradas que hice sobre el artículo Don’t Be Loopy: Re-Sampling and Simulation the SAS® Way.

Trata sobre lo siguiente:

Construir un cojunto de datos simples (dos vectores, x e y).
Hacer una regresión de y sobre x y capturar los residuos.
Crear 1000 vectores y' distintos añadiendo a $\hat{y}$ (la predicción de y) en el modelo anterior una reordenación de los residuos.
Crear los correspondientes 1000 modelos haciendo la regresión de cada $\hat{y}$ sobre x.
Obtener el histograma del coeficiente de la regresión.

Es un caso de bootstrap en el que no se muestrean directamente los valores iniciales sino los residuos del modelo.

El código en SAS es el siguiente:

 
/* 0: data creation */ 
data temp1;
	x=1; y=45; output;
	do x = 2 to 29;
		y = 3*x + 6*rannor(1234);
		output;
	end;
	x=30; y=45; output;
run;

%let regressors = x;
%let indata = temp1;

/* 1: perform the regression and get the predicted and residual values */ 
proc reg data= &INDATA;
	model y=®ressors;
	output out=out1 p=yhat r=res;
run;

/* 2: split the data: only the residuals will require URS */ 
data fit(keep=yhat ®RESSORS order) resid(keep=res);
	set out1;
	order+1;
run;

/* 3: this doesn't do any sampling - it copies the FIT data set repeatedly */ 
proc surveyselect data=fit out=outfit method=srs samprate=1 rep=1000;
run;

/* 4: this does the WR sampling of residuals for each replicate */ 
data outres2;
	do replicate = 1 to 1000;
		do order = 1 to numrecs;
			p = ceil(numrecs * ranuni(394747373));
			set resid nobs=numrecs point=p;
			output;
		end;
	end;
	stop;
run;

/* 5: then the randomized residuals are merged with the unrandomized records */ 
data prepped;
	merge outfit outres2;
	by replicate order;
	new_y=yhat+res;
run;

/* 6: the bootstrap process runs on each replicate */ 
proc reg data=prepped outest=est1(drop=_:);
	model new_y=®ressors;
	by replicate;
run;

/* 7: and the sampling distribution is aggregated */
proc univariate data=est1;
	var x;
	output out=final pctlpts=2.5, 97.5 pctlpre=ci;
run;

proc print;
run;

que corre en mi máquina en 37.37 segundos.

Ofrezco dos alternativas sustancialmente más sucintas en R. La primera es una reinterpretación literal en R del código anterior,

x <- 2:29
y <- 3 * x + 6 * rnorm( length(x) )
x <- c( 1, x, 30 )
y <- c( 45, y, 45 )

m0 <- lm( y ~ x )
yhat <- m0$fitted.values
res  <- m0$residuals

resultados <- replicate( 1000, lm( yhat + sample( res ) ~ x )$coefficients[["x"]] )
hist( resultados )

que corre en 2.46 segundos.

El segundo utiliza el paquete boot,

library( boot )

datos <- data.frame( y = yhat, x = x, res = res )

foo.boot <- function( datos, ind ){
 lm(y + res[ind] ~ x, data = datos)$coefficients[["x"]]
}

res <- boot( datos, foo.boot, R = 1000 )
plot( res )

que corre en 2.71 segundos.

El paquete boot es un candidato ideal para reimplementarlo utilizando paralelismo. Sin embargo, su mantenedor no parece estar por la labor. ¡Qué pena!

12 comentarios sobre “Don’t be loopy! (II)”

ffernandez 23 septiembre, 2011 10:52

¿Qué equipo tienes? A mi el primer código en SAS me tarda en ejecutar 1.46 segundos. El segundo código en R 2.34 segundos, el tercero 2.54… me dan ganas de revisar también «don’t by loopy (I)»…
ffernandez 23 septiembre, 2011 11:07

En el caso de (I) el código SAS tarda 58 segundos en mi equipo y 7 segundos el de R. Empate!!
datanalytics 23 septiembre, 2011 15:17

@ffernandez
Pues tengo un equipo no particularmente malo. He vuelto a correr el código y copio la parte más relevante:

NOTE: Interactivity disabled with BY processing.
NOTE: PROCEDURE REG used (Total process time):
real time 36.25 seconds
cpu time 3.67 seconds

NOTE: The data set WORK.EST1 has 1000 observations and 4 variables.

Mientras corre oigo sufrir a mi disco duro. ¿Será una cuestión de caché? Porque el tiempo de CPU es bajo.
Diego 23 septiembre, 2011 15:40

Con la laptop del curro , una lenovo thinkpad Intel Core i5 M520 @2.4 Gz , 3GB RAM el codigo de SAS corrio en real time 7.9 segundos (6.3 segundos CPU Time, la mayor parte debido a la regresión) y el de R 2.7 segundos. Saludos!!
Diego 23 septiembre, 2011 15:49

Para 10 mil vectores y 10 mil modelos, la diferencia se reduce, SAS aprox 74.7 segundos y la de R 40.22 segundos. Saludos!!!
datanalytics 23 septiembre, 2011 16:03

@Diego
Pues yo, con 10000 tengo que limpiar la pantalla del «log» a mano de vez en cuando porque se me llena.
ffernandez 23 septiembre, 2011 16:16

NOTE: Interactivity disabled with BY processing.
NOTE: PROCEDURE REG used (Total process time):
real time 1.56 seconds
cpu time 1.18 seconds

NOTE: The data set WORK.EST1 has 1000 observations and 4 variables.

En un Core i5 M560 @2.67 GHz 4GB RAM
ffernandez 23 septiembre, 2011 16:25

No se que ha cambiado en mi laptop de hace unos minutos para acá:
NOTE: Interactivity disabled with BY processing.
NOTE: PROCEDURE REG used (Total process time):
real time 0.04 seconds
cpu time 0.06 seconds

NOTE: The data set WORK.EST1 has 1000 observations and 4 variables.
ffernandez 23 septiembre, 2011 16:25

Ah, no, no ha sido mi laptop, he puesto un noprint en el reg… jejeje
datanalytics 26 septiembre, 2011 20:17

@ffernandez
¿Ha corrido en 0.06 segundos? ¿Cómo haces para que SAS no rasque el disco duro todo el tiempo? En mi equipo el problema no es la CPU sino el hecho de que SAS no para de acceder al disco.
ffernandez 27 septiembre, 2011 18:49

Pues conscientemente no he hecho nada… también es cierto que no he hecho yo la instalación, lo he probado en otro equipo con resultados similares…
Pingback: datanalytics » Dont be loopy! (III: jackknife y paralelismo)

Los comentarios están desabilitados.