HierRegressions
From FarmShare
(Difference between revisions)
Line 28: | Line 28: | ||
- | /farmshare/software/examples/R/hierregress/cheese.rda | + | == running on CPU == |
+ | |||
+ | save source code as rpudtestcpu.r | ||
+ | |||
+ | <source lang="r"> | ||
+ | library(bayesm) | ||
+ | |||
+ | data("cheese") | ||
+ | str("cheese") | ||
+ | |||
+ | retailer <- levels(cheese$RETAILER) | ||
+ | nreg <- length(retailer); nreg | ||
+ | |||
+ | regdata <- NULL | ||
+ | |||
+ | for (i in 1:nreg) { | ||
+ | filter <- cheese$RETAILER==retailer[i] | ||
+ | y <- log(cheese$VOLUME[filter]) | ||
+ | X <- cbind(1, # intercept placeholder | ||
+ | cheese$DISP[filter], | ||
+ | log(cheese$PRICE[filter])) | ||
+ | regdata[[i]] <- list(y=y, X=X) | ||
+ | } | ||
+ | |||
+ | Data <- list(regdata=regdata) | ||
+ | Mcmc <- list(R=2000) | ||
+ | |||
+ | system.time( | ||
+ | outcpu <- bayesm::rhierLinearModel( | ||
+ | Data=Data, | ||
+ | Mcmc=Mcmc)) | ||
+ | |||
+ | beta.3cpu <- mean(as.vector(outcpu$betadraw[, 3, 201:2000])) | ||
+ | beta.3cpu | ||
+ | </source> | ||
+ | |||
+ | run as: | ||
+ | |||
+ | <source lang="sh"> | ||
+ | cp -p /farmshare/software/examples/R/hierregress/cheese.rda . | ||
+ | R --no-save < rpudtestcpu.r | ||
+ | </source> | ||
+ | |||
+ | you should see output similar to: | ||
+ | |||
+ | <source lang="sh"> | ||
+ | bishopj@rye01:~$ R --no-save < rpudtestcpu.r | ||
+ | |||
+ | R version 2.15.2 (2012-10-26) -- "Trick or Treat" | ||
+ | Copyright (C) 2012 The R Foundation for Statistical Computing | ||
+ | ISBN 3-900051-07-0 | ||
+ | Platform: x86_64-pc-linux-gnu (64-bit) | ||
+ | |||
+ | R is free software and comes with ABSOLUTELY NO WARRANTY. | ||
+ | You are welcome to redistribute it under certain conditions. | ||
+ | Type 'license()' or 'licence()' for distribution details. | ||
+ | |||
+ | R is a collaborative project with many contributors. | ||
+ | Type 'contributors()' for more information and | ||
+ | 'citation()' on how to cite R or R packages in publications. | ||
+ | |||
+ | Type 'demo()' for some demos, 'help()' for on-line help, or | ||
+ | 'help.start()' for an HTML browser interface to help. | ||
+ | Type 'q()' to quit R. | ||
+ | |||
+ | [Previously saved workspace restored] | ||
+ | |||
+ | > library(bayesm) | ||
+ | > | ||
+ | > data("cheese") | ||
+ | > str("cheese") | ||
+ | chr "cheese" | ||
+ | > | ||
+ | > retailer <- levels(cheese$RETAILER) | ||
+ | > nreg <- length(retailer); nreg | ||
+ | [1] 88 | ||
+ | > | ||
+ | > regdata <- NULL | ||
+ | > | ||
+ | > for (i in 1:nreg) { | ||
+ | + filter <- cheese$RETAILER==retailer[i] | ||
+ | + y <- log(cheese$VOLUME[filter]) | ||
+ | + X <- cbind(1, # intercept placeholder | ||
+ | + cheese$DISP[filter], | ||
+ | + log(cheese$PRICE[filter])) | ||
+ | + regdata[[i]] <- list(y=y, X=X) | ||
+ | + } | ||
+ | > | ||
+ | > Data <- list(regdata=regdata) | ||
+ | > Mcmc <- list(R=2000) | ||
+ | > | ||
+ | > system.time( | ||
+ | + outcpu <- bayesm::rhierLinearModel( | ||
+ | + Data=Data, | ||
+ | + Mcmc=Mcmc)) | ||
+ | Z not specified -- putting in iota | ||
+ | |||
+ | Starting Gibbs Sampler for Linear Hierarchical Model | ||
+ | 88 Regressions | ||
+ | 1 Variables in Z (if 1, then only intercept) | ||
+ | |||
+ | Prior Parms: | ||
+ | Deltabar | ||
+ | [,1] [,2] [,3] | ||
+ | [1,] 0 0 0 | ||
+ | A | ||
+ | [,1] | ||
+ | [1,] 0.01 | ||
+ | nu.e (d.f. parm for regression error variances)= 3 | ||
+ | Vbeta ~ IW(nu,V) | ||
+ | nu = 6 | ||
+ | V | ||
+ | [,1] [,2] [,3] | ||
+ | [1,] 6 0 0 | ||
+ | [2,] 0 6 0 | ||
+ | [3,] 0 0 6 | ||
+ | |||
+ | MCMC parms: | ||
+ | R= 2000 keep= 1 | ||
+ | |||
+ | MCMC Iteration (est time to end - min) | ||
+ | 100 ( 0.7 ) | ||
+ | 200 ( 0.7 ) | ||
+ | 300 ( 0.6 ) | ||
+ | 400 ( 0.6 ) | ||
+ | 500 ( 0.5 ) | ||
+ | 600 ( 0.5 ) | ||
+ | 700 ( 0.5 ) | ||
+ | 800 ( 0.4 ) | ||
+ | 900 ( 0.4 ) | ||
+ | 1000 ( 0.4 ) | ||
+ | 1100 ( 0.3 ) | ||
+ | 1200 ( 0.3 ) | ||
+ | 1300 ( 0.2 ) | ||
+ | 1400 ( 0.2 ) | ||
+ | 1500 ( 0.2 ) | ||
+ | 1600 ( 0.1 ) | ||
+ | 1700 ( 0.1 ) | ||
+ | 1800 ( 0.1 ) | ||
+ | 1900 ( 0 ) | ||
+ | 2000 ( 0 ) | ||
+ | Total Time Elapsed: 0.71 | ||
+ | user system elapsed | ||
+ | 42.288 0.036 42.448 | ||
+ | > | ||
+ | > beta.3cpu <- mean(as.vector(outcpu$betadraw[, 3, 201:2000])) | ||
+ | > beta.3cpu | ||
+ | [1] -2.146799 | ||
+ | > | ||
+ | </source> |
Revision as of 10:43, 8 September 2013
Hierarchical Linear Regressions
http://www.r-tutor.com/gpu-computing/rbayes/rhierlmc
http://www.r-tutor.com/content/download
installing rpud
module load cuda mkdir /tmp/$USER cd /tmp/$USER wget http://www.r-tutor.com/sites/default/files/rpud/rpud_0.3.3.tar.gz mkdir ~/lib/R export R_LIBS_USER="${HOME}/lib/R" NVCC_OPT="-arch compute_20" R CMD INSTALL --library=$R_LIBS rpud_0.3.3.tar.gz
module load cuda export R_LIBS_USER="${HOME}/lib/R" wget http://www.r-tutor.com/sites/default/files/rpud/rpudplus_0.3.3.tar.gz cd rpudplus_0.3.3/ ./install.sh
running on CPU
save source code as rpudtestcpu.r
library(bayesm) data("cheese") str("cheese") retailer <- levels(cheese$RETAILER) nreg <- length(retailer); nreg regdata <- NULL for (i in 1:nreg) { filter <- cheese$RETAILER==retailer[i] y <- log(cheese$VOLUME[filter]) X <- cbind(1, # intercept placeholder cheese$DISP[filter], log(cheese$PRICE[filter])) regdata[[i]] <- list(y=y, X=X) } Data <- list(regdata=regdata) Mcmc <- list(R=2000) system.time( outcpu <- bayesm::rhierLinearModel( Data=Data, Mcmc=Mcmc)) beta.3cpu <- mean(as.vector(outcpu$betadraw[, 3, 201:2000])) beta.3cpu
run as:
cp -p /farmshare/software/examples/R/hierregress/cheese.rda . R --no-save < rpudtestcpu.r
you should see output similar to:
bishopj@rye01:~$ R --no-save < rpudtestcpu.r R version 2.15.2 (2012-10-26) -- "Trick or Treat" Copyright (C) 2012 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. [Previously saved workspace restored] > library(bayesm) > > data("cheese") > str("cheese") chr "cheese" > > retailer <- levels(cheese$RETAILER) > nreg <- length(retailer); nreg [1] 88 > > regdata <- NULL > > for (i in 1:nreg) { + filter <- cheese$RETAILER==retailer[i] + y <- log(cheese$VOLUME[filter]) + X <- cbind(1, # intercept placeholder + cheese$DISP[filter], + log(cheese$PRICE[filter])) + regdata[[i]] <- list(y=y, X=X) + } > > Data <- list(regdata=regdata) > Mcmc <- list(R=2000) > > system.time( + outcpu <- bayesm::rhierLinearModel( + Data=Data, + Mcmc=Mcmc)) Z not specified -- putting in iota Starting Gibbs Sampler for Linear Hierarchical Model 88 Regressions 1 Variables in Z (if 1, then only intercept) Prior Parms: Deltabar [,1] [,2] [,3] [1,] 0 0 0 A [,1] [1,] 0.01 nu.e (d.f. parm for regression error variances)= 3 Vbeta ~ IW(nu,V) nu = 6 V [,1] [,2] [,3] [1,] 6 0 0 [2,] 0 6 0 [3,] 0 0 6 MCMC parms: R= 2000 keep= 1 MCMC Iteration (est time to end - min) 100 ( 0.7 ) 200 ( 0.7 ) 300 ( 0.6 ) 400 ( 0.6 ) 500 ( 0.5 ) 600 ( 0.5 ) 700 ( 0.5 ) 800 ( 0.4 ) 900 ( 0.4 ) 1000 ( 0.4 ) 1100 ( 0.3 ) 1200 ( 0.3 ) 1300 ( 0.2 ) 1400 ( 0.2 ) 1500 ( 0.2 ) 1600 ( 0.1 ) 1700 ( 0.1 ) 1800 ( 0.1 ) 1900 ( 0 ) 2000 ( 0 ) Total Time Elapsed: 0.71 user system elapsed 42.288 0.036 42.448 > > beta.3cpu <- mean(as.vector(outcpu$betadraw[, 3, 201:2000])) > beta.3cpu [1] -2.146799 >