HierRegressions

From FarmShare

(Difference between revisions)
Jump to: navigation, search
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
>
Retrieved from "https://web.stanford.edu/group/farmshare/cgi-bin/wiki/index.php/HierRegressions"
Personal tools
Toolbox
LANGUAGES