HierRegressions
From FarmShare
(Difference between revisions)
Line 176: | Line 176: | ||
> beta.3cpu | > beta.3cpu | ||
[1] -2.146799 | [1] -2.146799 | ||
+ | > | ||
+ | </source> | ||
+ | |||
+ | == running on GPU == | ||
+ | |||
+ | |||
+ | <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) | ||
+ | |||
+ | library(rpud) # load rpudplus | ||
+ | |||
+ | system.time( | ||
+ | outgpu <- rpud::rhierLinearModel( | ||
+ | Data=Data, | ||
+ | Mcmc=Mcmc, | ||
+ | output="bayesm")) | ||
+ | |||
+ | beta.3 <- mean(as.vector(outgpu$betadraw[, 3, 201:2000])) | ||
+ | beta.3 | ||
+ | </source> | ||
+ | |||
+ | <source lang="sh"> | ||
+ | bishopj@rye01:~$ module load cuda | ||
+ | bishopj@rye01:~$ export R_LIBS_USER=$HOME/lib/R | ||
+ | bishopj@rye01:~$ R --no-save < rpudtestgpu.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) | ||
+ | > | ||
+ | > library(rpud) # load rpudplus | ||
+ | Rpudplus 0.3.3 | ||
+ | http://www.r-tutor.com | ||
+ | Copyright (C) 2010-2013 Chi Yau. All Rights Reserved. | ||
+ | Rpudplus is free for academic use only. There is absolutely NO warranty. | ||
+ | |||
+ | Attaching package: 'rpud' | ||
+ | |||
+ | The following object(s) are masked from 'package:bayesm': | ||
+ | |||
+ | rhierLinearModel, rmultireg | ||
+ | |||
+ | > | ||
+ | > system.time( | ||
+ | + outgpu <- rpud::rhierLinearModel( | ||
+ | + Data=Data, | ||
+ | + Mcmc=Mcmc, | ||
+ | + output="bayesm")) | ||
+ | |||
+ | 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. 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 | ||
+ | |||
+ | .................... | ||
+ | user system elapsed | ||
+ | 0.192 0.156 0.350 | ||
+ | > | ||
+ | > beta.3 <- mean(as.vector(outgpu$betadraw[, 3, 201:2000])) | ||
+ | > beta.3 | ||
+ | [1] -2.145573 | ||
> | > | ||
</source> | </source> |
Revision as of 10:46, 8 September 2013
Contents |
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 >
running on GPU
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) library(rpud) # load rpudplus system.time( outgpu <- rpud::rhierLinearModel( Data=Data, Mcmc=Mcmc, output="bayesm")) beta.3 <- mean(as.vector(outgpu$betadraw[, 3, 201:2000])) beta.3
bishopj@rye01:~$ module load cuda bishopj@rye01:~$ export R_LIBS_USER=$HOME/lib/R bishopj@rye01:~$ R --no-save < rpudtestgpu.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) > > library(rpud) # load rpudplus Rpudplus 0.3.3 http://www.r-tutor.com Copyright (C) 2010-2013 Chi Yau. All Rights Reserved. Rpudplus is free for academic use only. There is absolutely NO warranty. Attaching package: 'rpud' The following object(s) are masked from 'package:bayesm': rhierLinearModel, rmultireg > > system.time( + outgpu <- rpud::rhierLinearModel( + Data=Data, + Mcmc=Mcmc, + output="bayesm")) 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. 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 .................... user system elapsed 0.192 0.156 0.350 > > beta.3 <- mean(as.vector(outgpu$betadraw[, 3, 201:2000])) > beta.3 [1] -2.145573 >