Multi-dimensional data and parameters: Is this the most efficient way to code this problem?
227 views
Skip to first unread message
HookahBoy
Nov 24, 2012, 10:41:25 AM11/24/12
to stan-...@googlegroups.com, Gyan Sinha
Dear All:
I am still a bit fuzzy about arrays, vectors, and matrices in Stan. In particular, some of the issues relating to what final form results from an operation on any two of these basic types and also how Stan is structured to prefer one form over the other in terms of efficiency.
I am attaching an example that is working for me (extension of the 8.schools example) but I am not sure I have coded it in the most efficient manner possible, so here goes:
#####################################
require(rstan)
require(plyr)
require(doMC)
require(coda)
require(ggmcmc)
options(cores=4)
registerDoMC()
norm_code <- '
data {
int<lower=0> N;
int<lower=0> D;
row_vector[D] Y[N];
row_vector<lower=0>[D] S[N];
}
parameters {
real mu[D];
real<lower=0> tau[D];
row_vector[D] eta[N];
}
transformed parameters {
row_vector[D] theta[N];
for (n in 1:N) {
for (d in 1:D)
theta[n,d]<- mu[d] + tau[d] * eta[n,d];
}
}
model {
for (n in 1:N)
for (d in 1:D) {
eta[n,d] ~ normal(0, 1);
Y[n,d] ~ normal(theta[n,d], S[n,d]);
}
}
'
## Remove NAs from data and choose shelf
### which.ones <- !is.na(info.param.p[,'age']) & info.param.p[,'shelf'] == 'AHM'
### Y <- data.matrix(info.param.p[which.ones,-c(1:3)])
### S <- data.matrix(info.se.p[which.ones,-c(1:3)])
Y <- read.csv('Y.csv')
S <- read.csv('S.csv')
N <- dim(Y)[1]
D <- dim(Y)[2]
norm_dat <- list(Y=Y,S=S,N=N,D=D)
initial_multi_lvl_fit <- stan(model_code = norm_code, data = norm_dat, iter = 1,
chains = 1)
num_chains <- 4
num_iter <- 5000
#prep for parallel case
fun <- function(){
fit <- stan(
fit = initial_multi_lvl_fit
, data = norm_dat
, iter = num_iter,
, chains = 1
, verbose = F
)
return(fit)
}
#parallel case
start <- proc.time()[3]
parallel_fit <- foreach(i = 1:num_chains) %dopar% fun()
do.plots <- FALSE
s.try <- extract(parallel_fit[[1]])
params.to.exclude <- -c((D*2 + 1):(D*2 + N*D),dim(s.try)[3])
s <- array(dim=c(dim(s.try)[1],num_chains,dim(s.try)[3]),
dimnames=list('iterations'=NULL,'chains'=paste('chain:',1:num_chains,sep=''),
'parameters'=dimnames(s.try)$parameters))
if (do.plots) {
for (i in 1:num_chains) {
s[,i,] <- extract(parallel_fit[[i]])
}
s.list <- do.call(mcmc.list,alply(s[,,params.to.exclude],2,mcmc))
S <- suppressWarnings(ggs(s.list))
ggmcmc(S, file = 'multi.lvl.V2.pdf')
}
theta.smp <- s.try[,,c(-(1:10),params.to.exclude)]
theta.df <- matrix(0,nrow=N,ncol=D)
j <- 1
for (i in 1:D) {
cat('j:j+N ',j,j+N - 1,'\n')
theta.df[,i] <- apply(theta.smp[,j:(j+N - 1)],2,median)
j <- j + N
}
rownames(theta.df) <- info.param.p[rownames(Y),'pool.id']
colnames(theta.df) <- colnames(Y)
theta.df <- data.frame(pool.id = rownames(theta.df),theta.df)
rownames(theta.df) <- NULL
elapsed.time <- proc.time()[3] - start
print(elapsed.time)
################ sessionInfo() ########
R version 2.15.0 (2012-03-30)
Platform: x86_64-pc-linux-gnu (64-bit)
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=C LC_COLLATE=C LC_MONETARY=C LC_MESSAGES=C
[7] LC_PAPER=C LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=C LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] ggmcmc_0.2 ggplot2_0.9.2.1 reshape_0.8.4 coda_0.14-7 doMC_1.2.5 multicore_0.1-7 foreach_1.3.5 codetools_0.2-8
[9] iterators_1.0.5 plyr_1.7.1 rstan_1.0.3 RcppEigen_0.3.1 Matrix_1.0-6 lattice_0.20-6 inline_0.3.10 Rcpp_0.9.10
loaded via a namespace (and not attached):
[1] MASS_7.3-18 RColorBrewer_1.0-5 colorspace_1.1-1 compiler_2.15.0 dichromat_1.2-4 digest_0.5.2 grid_2.15.0
[8] gtable_0.1.1 labeling_0.1 memoise_0.1 munsell_0.3 proto_0.3-9.2 reshape2_1.2.1 scales_0.2.2
[15] stringr_0.6 tools_2.15.0
I am still a bit fuzzy about arrays, vectors, and matrices in Stan. In particular, some of the issues relating to what final form results from an operation on any two of these basic types and also how Stan is structured to prefer one form over the other in terms of efficiency.
I am attaching an example that is working for me (extension of the 8.schools example) but I am not sure I have coded it in the most efficient manner possible, so here goes:
#####################################
require(rstan)
require(plyr)
require(doMC)
require(coda)
require(ggmcmc)
options(cores=4)
registerDoMC()
norm_code <- '
data {
int<lower=0> N;
int<lower=0> D;
row_vector[D] Y[N];
row_vector<lower=0>[D] S[N];
}
parameters {
real mu[D];
real<lower=0> tau[D];
row_vector[D] eta[N];
}
transformed parameters {
row_vector[D] theta[N];
for (n in 1:N) {
for (d in 1:D)
theta[n,d]<- mu[d] + tau[d] * eta[n,d];
}
}
model {
for (n in 1:N)
for (d in 1:D) {
eta[n,d] ~ normal(0, 1);
Y[n,d] ~ normal(theta[n,d], S[n,d]);
}
}
'
## Remove NAs from data and choose shelf
### which.ones <- !is.na(info.param.p[,'age']) & info.param.p[,'shelf'] == 'AHM'
### Y <- data.matrix(info.param.p[which.ones,-c(1:3)])
### S <- data.matrix(info.se.p[which.ones,-c(1:3)])
Y <- read.csv('Y.csv')
S <- read.csv('S.csv')
N <- dim(Y)[1]
D <- dim(Y)[2]
norm_dat <- list(Y=Y,S=S,N=N,D=D)
initial_multi_lvl_fit <- stan(model_code = norm_code, data = norm_dat, iter = 1,
chains = 1)
num_chains <- 4
num_iter <- 5000
#prep for parallel case
fun <- function(){
fit <- stan(
fit = initial_multi_lvl_fit
, data = norm_dat
, iter = num_iter,
, chains = 1
, verbose = F
)
return(fit)
}
#parallel case
start <- proc.time()[3]
parallel_fit <- foreach(i = 1:num_chains) %dopar% fun()
do.plots <- FALSE
s.try <- extract(parallel_fit[[1]])
params.to.exclude <- -c((D*2 + 1):(D*2 + N*D),dim(s.try)[3])
s <- array(dim=c(dim(s.try)[1],num_chains,dim(s.try)[3]),
dimnames=list('iterations'=NULL,'chains'=paste('chain:',1:num_chains,sep=''),
'parameters'=dimnames(s.try)$parameters))
if (do.plots) {
for (i in 1:num_chains) {
s[,i,] <- extract(parallel_fit[[i]])
}
s.list <- do.call(mcmc.list,alply(s[,,params.to.exclude],2,mcmc))
S <- suppressWarnings(ggs(s.list))
ggmcmc(S, file = 'multi.lvl.V2.pdf')
}
theta.smp <- s.try[,,c(-(1:10),params.to.exclude)]
theta.df <- matrix(0,nrow=N,ncol=D)
j <- 1
for (i in 1:D) {
cat('j:j+N ',j,j+N - 1,'\n')
theta.df[,i] <- apply(theta.smp[,j:(j+N - 1)],2,median)
j <- j + N
}
rownames(theta.df) <- info.param.p[rownames(Y),'pool.id']
colnames(theta.df) <- colnames(Y)
theta.df <- data.frame(pool.id = rownames(theta.df),theta.df)
rownames(theta.df) <- NULL
elapsed.time <- proc.time()[3] - start
print(elapsed.time)
################ sessionInfo() ########
R version 2.15.0 (2012-03-30)
Platform: x86_64-pc-linux-gnu (64-bit)
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=C LC_COLLATE=C LC_MONETARY=C LC_MESSAGES=C
[7] LC_PAPER=C LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=C LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] ggmcmc_0.2 ggplot2_0.9.2.1 reshape_0.8.4 coda_0.14-7 doMC_1.2.5 multicore_0.1-7 foreach_1.3.5 codetools_0.2-8
[9] iterators_1.0.5 plyr_1.7.1 rstan_1.0.3 RcppEigen_0.3.1 Matrix_1.0-6 lattice_0.20-6 inline_0.3.10 Rcpp_0.9.10
loaded via a namespace (and not attached):
[1] MASS_7.3-18 RColorBrewer_1.0-5 colorspace_1.1-1 compiler_2.15.0 dichromat_1.2-4 digest_0.5.2 grid_2.15.0
[8] gtable_0.1.1 labeling_0.1 memoise_0.1 munsell_0.3 proto_0.3-9.2 reshape2_1.2.1 scales_0.2.2
[15] stringr_0.6 tools_2.15.0
Ben Goodrich
Nov 24, 2012, 1:26:35 PM11/24/12
to stan-...@googlegroups.com, Gyan Sinha
On Saturday, November 24, 2012 10:41:25 AM UTC-5, HookahBoy wrote:
I am still a bit fuzzy about arrays, vectors, and matrices in Stan. In particular, some of the issues relating to what final form results from an operation on any two of these basic types and also how Stan is structured to prefer one form over the other in terms of efficiency.
I would say the differences in efficiency are a minor consideration; the main thing is whether the .stan file is easy to understand. To those ends, I think three principles can be established:
1) Unless you have a compelling reason not to, use matrices because those are handled by Stan in pretty much the same way as in R.
2) If you use the "collection of vectors" approach in Stan, use the transformed data {} block to print() parts of your data to make sure they are the way you intended.
3) As of the not-yet-released RStan 1.0.4, if you use the "collection of vectors" approach in Stan, pass the data as a list of objects, which (in principle) will either do the right thing or give you an informative message telling you how to do the right thing.
All that said, I think you managed to write a valid .stan file in this case.
I am attaching an example that is working for me (extension of the 8.schools example) but I am not sure I have coded it in the most efficient manner possible, so here goes:
No matter how you do the input / output of data, you would want to vectorize the sampling statements as much as possible. Here is how I would have done it:
data {
int<lower=0> N;
int<lower=0> D;
matrix[N,D] Y;
matrix[N,D] S;
}
parameters {
vector[D] mu;
vector<lower=0>[D] tau;
matrix[N,D] eta;
}
transformed parameters {
matrix[N,D] theta;
for (n in 1:N) {
for (d in 1:D)
theta[n,d] <- mu[d] + tau[d] * eta[n,d];
}
}
model {
// this will work in Stan >= 1.0.4:
// eta ~ normal(0.0, 1.0);
// Y ~ normal(theta, S);
for (d in 1:D) {
col(eta,d) ~ normal(0, 1);
col(Y,d) ~ normal(col(theta,d), col(S,d));
}
}
Your R code to process the results seems more involved than is necessary. I think you can get away with
parallel_fit <- sflist2stanfit( foreach(i = 1:num_chains) %dopar% fun() )
s <- as.array(parallel_fit)
included_parameters <- c("something", "else")
s.list <- mcmc.list(apply(s[,,included_parameters], 2, mcmc))
and go from there.
Ben
HookahBoy
Nov 24, 2012, 7:54:20 PM11/24/12
to stan-...@googlegroups.com, Gyan Sinha
Ben:
On Saturday, November 24, 2012 10:41:25 AM UTC-5, HookahBoy wrote:
Thanks for your response and looking forward to the release of Stan/RStan 1.0.4. Can you give us a sense of when they might be out?
Best Regards
On Saturday, November 24, 2012 10:41:25 AM UTC-5, HookahBoy wrote:
Ben Goodrich
Nov 24, 2012, 10:43:22 PM11/24/12
to stan-...@googlegroups.com, Gyan Sinha
December
Bob Carpenter
Nov 25, 2012, 3:25:42 PM11/25/12
to stan-...@googlegroups.com
On 11/24/12 10:41 AM, HookahBoy wrote:
> I am still a bit fuzzy about arrays, vectors, and matrices in Stan. In particular, some of the issues relating to what
> final form results from an operation on any two of these basic types and also how Stan is structured to prefer one form
> over the other in terms of efficiency.
you do, it's negligible compared to the gradient calculations.
Even if you define things by component in a loop, the loops are
compiled.
Here's how to infer the type of an expression:
literal 1 integer
literal 1.0 real
variable x given by declaration (no hiding, so only one)
x[i1,...,iN] knock N indices off of type of x
(knocking one index off of matrix gives row vector,
one index off of vector or row_vector gives real)
f(e1,...,eN) figure out types of subexpressions
e1,...,eN, look up f in the appendix
and read off the result type
(int types may be promoted to real if necessary)
In looking up functions, binary operators like * are
defined as operator*() in the documentation and index.
Matrices declared as:
matrix[M,N] x[I,J];
are indexed as x[i,j,m,n] -- the array indices come first.
Also, all of the matrix operations are defined to return
the lowest inferrable type; for example,
row_vector * vector
produces a real, not a matrix.
I'll write a section to this effect in the manual so I
can just point to it in the future.
- Bob
Bob Carpenter
Nov 25, 2012, 5:43:54 PM11/25/12
to stan-...@googlegroups.com
Stan/RStan 1.0.4 is pretty much ready to go
other than some test issues. We may not have
time to do the release for a week or two given
our various travel schedules.
Here's what's in it so far (though I may have missed
something and the jury's still out on changing categorical()
support from 0:(N-1) to 1:N and on removing non-symmetric
eigenvalues() functions):
v1.0.4
======================================================================
-- Libraries Updated
* Eigen version 3.1.1 replaced with version 3.1.2
* Boost version 1.51.0 replaced with version 1.52.0
-- Bug Fixes
* fixed bug introduced in 1.0.3 that hid line and column number
bug reporting
* added checks that data dimensions match as well as sizes
* removed non-symmetric versions of eigenvalues() and eigenvectors()
-- Improvements
* generating double-based log-prob function (in addition to
stan::agrad::var-based)
* testing identifiers are not reserved words in C++/Stan
* trapping/reporting locations of errors in data and init reads
* templated out matrix libs to reduce code duplication
* generalized multiply_lower_tri_self_transpose() to non-square
matrices
* vectorized auto-dif for tcrossprod() and crossprod()
-- New Functions
* special functions: log_inv_logit(), log1m_inv_logit()
* matrix special functions: cumulative_sum()
* probability functions: poisson_log_log() for log-rate
parameterized Poisson
* matrix functions: block(), diag_pre_multiply(), diag_post_multiply()
- Bob
other than some test issues. We may not have
time to do the release for a week or two given
our various travel schedules.
Here's what's in it so far (though I may have missed
something and the jury's still out on changing categorical()
support from 0:(N-1) to 1:N and on removing non-symmetric
eigenvalues() functions):
v1.0.4
======================================================================
-- Libraries Updated
* Eigen version 3.1.1 replaced with version 3.1.2
* Boost version 1.51.0 replaced with version 1.52.0
-- Bug Fixes
* fixed bug introduced in 1.0.3 that hid line and column number
bug reporting
* added checks that data dimensions match as well as sizes
* removed non-symmetric versions of eigenvalues() and eigenvectors()
-- Improvements
* generating double-based log-prob function (in addition to
stan::agrad::var-based)
* testing identifiers are not reserved words in C++/Stan
* trapping/reporting locations of errors in data and init reads
* templated out matrix libs to reduce code duplication
* generalized multiply_lower_tri_self_transpose() to non-square
matrices
* vectorized auto-dif for tcrossprod() and crossprod()
-- New Functions
* special functions: log_inv_logit(), log1m_inv_logit()
* matrix special functions: cumulative_sum()
* probability functions: poisson_log_log() for log-rate
parameterized Poisson
* matrix functions: block(), diag_pre_multiply(), diag_post_multiply()
- Bob
Tim Triche, Jr.
Nov 26, 2012, 1:54:37 PM11/26/12
to stan-...@googlegroups.com
This looks awesome. Is it relatively safe to use what's in the Git repo on googlecode?
Bob Carpenter
Nov 27, 2012, 2:34:14 AM11/27/12
to stan-...@googlegroups.com
On 11/26/12 1:54 PM, Tim Triche, Jr. wrote:
> This looks awesome. Is it relatively safe to use what's in the Git repo on googlecode?
more stable :-)
We aren't super picky about making sure every commit
keeps absolutely everything working. If all the unit
and model tests pass, things should be OK.
Right now, the unit tests are messed up due to
some generalizations using C++11 features for
testing vectorized model gradients -- as soon as that's
sorted, we're doing the 1.0.4 release.
In the meantime, I've added (but not pushed)
definitions for boolean operators (!, &&, ||) and
comparison operators (<, >, <=, >=, ==, !=), returning
int types. It's mostly tested and should be ready for
1.0.4. I may even have time for if-else statements
and general while loops.
- Bob
>
>
> On Sun, Nov 25, 2012 at 2:43 PM, Bob Carpenter <ca...@alias-i.com <mailto:ca...@alias-i.com>> wrote:
>
> Stan/RStan 1.0.4 is pretty much ready to go
> other than some test issues. We may not have
> time to do the release for a week or two given
> our various travel schedules.
>
> Here's what's in it so far (though I may have missed
> something and the jury's still out on changing categorical()
> support from 0:(N-1) to 1:N and on removing non-symmetric
> eigenvalues() functions):
>
> v1.0.4
> -- Libraries Updated
> * Eigen version 3.1.1 replaced with version 3.1.2
> * Boost version 1.51.0 replaced with version 1.52.0
> -- Bug Fixes
> * fixed bug introduced in 1.0.3 that hid line and column number
> bug reporting
> * added checks that data dimensions match as well as sizes
> * removed non-symmetric versions of eigenvalues() and eigenvectors()
> -- Improvements
> * generating double-based log-prob function (in addition to
> stan::agrad::var-based)
> * testing identifiers are not reserved words in C++/Stan
> * trapping/reporting locations of errors in data and init reads
> * templated out matrix libs to reduce code duplication
> * generalized multiply_lower_tri_self___transpose() to non-square
> * Eigen version 3.1.1 replaced with version 3.1.2
> * Boost version 1.51.0 replaced with version 1.52.0
> -- Bug Fixes
> * fixed bug introduced in 1.0.3 that hid line and column number
> bug reporting
> * added checks that data dimensions match as well as sizes
> * removed non-symmetric versions of eigenvalues() and eigenvectors()
> -- Improvements
> * generating double-based log-prob function (in addition to
> stan::agrad::var-based)
> * testing identifiers are not reserved words in C++/Stan
> * trapping/reporting locations of errors in data and init reads
> * templated out matrix libs to reduce code duplication
> matrices
> * vectorized auto-dif for tcrossprod() and crossprod()
> -- New Functions
> * special functions: log_inv_logit(), log1m_inv_logit()
> * matrix special functions: cumulative_sum()
> * probability functions: poisson_log_log() for log-rate
> parameterized Poisson
> * matrix functions: block(), diag_pre_multiply(), diag_post_multiply()
>
> - Bob
>
> On 11/24/12 10:43 PM, Ben Goodrich wrote:
>
> December
>
> On Saturday, November 24, 2012 7:54:20 PM UTC-5, HookahBoy wrote:
>
> Ben:
>
> Thanks for your response and looking forward to the release of Stan/RStan 1.0.4. Can you give us a sense of
> when
> they might be out?
>
>
> --
>
>
>
>
>
> --
> /A model is a lie that helps you see the truth./
> * vectorized auto-dif for tcrossprod() and crossprod()
> -- New Functions
> * special functions: log_inv_logit(), log1m_inv_logit()
> * matrix special functions: cumulative_sum()
> * probability functions: poisson_log_log() for log-rate
> parameterized Poisson
> * matrix functions: block(), diag_pre_multiply(), diag_post_multiply()
>
> - Bob
>
> On 11/24/12 10:43 PM, Ben Goodrich wrote:
>
> December
>
> On Saturday, November 24, 2012 7:54:20 PM UTC-5, HookahBoy wrote:
>
> Ben:
>
> Thanks for your response and looking forward to the release of Stan/RStan 1.0.4. Can you give us a sense of
> when
> they might be out?
>
>
> --
>
>
>
>
>
> --
> /
> /
> Howard Skipper <http://cancerres.aacrjournals.org/content/31/9/1173.full.pdf>
>
> --
>
>
Reply all
Reply to author
Forward
0 new messages