feature request for rstan

[Sebastian Weber]

Hi!

 

I am regurlary handling models where I need to know during initialization the number of groups which are fitted in my hierachical model. This is needed to seed the random effects defined in the model for the individual units. For that to work out the init functions needs to figure some things out about the data which I solve right now like this

 

attach(stan_data)

stan("mymodel.stan", data=stan_data, init=myinitFunction)

detach(stan_data)

 

Assume now that stan_data$J <- 20 is the number of individual units over which random effects are defined. Then the init function should look like

 

myinitFunction <- function() { list(randomEffect=rnorm(J,0,0.5)) }

 

My trick with attach works, but I would certainly consider this as a hack. So if someone else finds this useful, I would like to request that the data is made available to the myinitFunction when this function is called.

 

I hope this is not to confused - I am happy to work out a minimal example which even works.

 

Best,

 

Sebastian

[Bob Carpenter]

We're about to undertake a major overhaul on the RStan interface,
so now is the time to get feature requests in!

Sebastian's "trick" is exploiting Stan's dynamic scope to
read a global variable in the body of a function. Perhaps not
the safest programming style, but then I'm not a huge fan of
this kind of scoping.

But R is nothing if not flexible. Maybe you'd be happier with
a functional programming solution --- take the data as an argument
to a "factory" function that creates the init function. For example,

init_fun_factory <- function(stan_data) {
return(function() { list(randomEffect=rnorm(stan_data$J,0,0.5)); }
}

and then call with

stan("mymodel.stan", data=stan_data, init=init_fun_factory(stan_data));

Here's a simple example of how the idea works outside the context
of Stan:

> factory <- function(data) return(function() { return(data$alpha); })

> data <- list(alpha=10)

> alpha
Error: object 'alpha' not found

> init_fun <- factory(data)

> init_fun()
[1] 10

You can also just use factory(data)() to get 10.

Was there something you were suggesting we add to rstan?

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Sebastian Weber]

Hi again!

 

Just learned once more how nasty R can be: One can change the calling enviroment of a function, such that this code snippet

 

    if("init" %in% names(stan_args)) {
        if(is.function(stan_args$init))
            environment(stan_args$init) <- list2env(stan_data)
    }
    fit <- do.call(stan_run, stan_args)

Does what I want, i.e. if init is a function, then I make the data the encapuslating environment for init. Maybe this could become the default for rstan?

 

Best,

Sebastian

[Bob Carpenter]

My own preference is not to have complex defaults and not under
any circumstances to play with default environments for user-supplied
functions. The former so everything's easy to doc and test and understand
(though perhaps not so easy to use in a single step) and the latter
so that the behavior is predictible to R users who might themselves
be monkeying with environments as you did.

But I'm not the only one with my hands on the steering wheel of the
Stan bus. And I'm more concerned that the lower-level functionality is
implemented cleanly than the high-level calls (for instance, "stan()" calls
compilers, data readers, and then one or more chain runners --- it's
the thing stan() calls that I'm most concerned with).

That %in% function looks really useful!

- Bob

[Ben Goodrich]

On Wednesday, July 9, 2014 12:46:32 PM UTC-4, Bob Carpenter wrote:

But I'm not the only one with my hands on the steering wheel of the
Stan bus.  And I'm more concerned that the lower-level functionality is
implemented cleanly than the high-level calls (for instance, "stan()" calls
compilers, data readers, and then one or more chain runners --- it's
the thing stan() calls that I'm most concerned with).

I'm inclined to agree that we not mess with this right now. Under my proposal for 3.0, there could be an abstract base class for the stanfit object and then we could create an child class on the fly that would have additional fields for the symbols declared in the data block of the .stan program that get populated (by reference) with the elements of the passed data list in R. Then, it wouldn't be too difficult for a function that generates initial values to refer internally to "this", as in

this$get_J()

Ben

[Krzysztof Sakrejda]

On Wednesday, July 9, 2014 12:38:32 PM UTC-4, Bob Carpenter wrote:

We're about to undertake a major overhaul on the RStan interface,
so now is the time to get feature requests in!

Since you asked...

Feature request: I would like to be able to ask Stan to sample p(x|y,data) and use other code to sample p(y|x,data) to get samples from p(x,y|data). I think what is needed to make that happen is:
   1) make it easy to read/write data and read/write parameter values at any time from R without otherwise modifying the state/adaptation of the sampler.
   2) make it easy to read/write the sampler/optimizer tuning parameters from R without modifying the rest of the object, redoing adaptation, etc...

Argument:
- Stan can do a lot but it can't do everything. 
- R is frustrating at times but it's a good glue language for statistical analysis because everyone feels obliged to implement interfaces there.
- I acknowledge that it is sometimes necessary to re-do the adaptation but of the many ways you give users to shoot themselves in the foot this would not be the worst. 
- Right now when I do an analysis I have to ask myself whether it's easier better to write my own code OR use Stan, I would like to be able to make that an "AND". 
- This would provide a path for people outside of the Stan team to show that certain algorithms/functionality are relevant for Stan without writing Stan-compatible C++ code.
- Working through R would be relatively slow but in my experience most of the wait time is due to autocorrelation caused by the model/data rather than the speed of the algorithm.

Not my first time pitching this one but thought I'd give it another try...

Krzysztof

 

[Bob Carpenter]

Let me rewrite your request following Andrew's (and hence
Stan's) notation.

You have parameters theta = (alpha,beta) and data y and
you want to sample from a posterior p(alpha,beta | y)
by doing something Gibbs-like in sampling p(alpha | beta, y)
in Stan and p(beta | alpha, y) in some other program.

Fair enough --- we'd love to make that easy to do. In fact,
it was part of our rejected DARPA probabilistic programming
proposal.

I think we'll be adding just about everything you'll need in
the next version (more below).

WARNING: I think this is what you meant below about adaptation,
but for everyone else's sake: because beta now looks like "data" to
Stan (even though it's a parameter of the bigger system), every time beta
changes, it changes the posterior from which Stan is trying to sample.
This is an issue because our adaptation assumes that we're
only trying to adapt to a single posterior.

In R, we are going to be breaking the stan() command down and exposing
some more primitive operations:

* compiling a model: translates model from Stan to C++ and compiles
a dynamically linkable library, but doesn't instantiate any instances
of the model

* create a model object with data: populate a compiled model with data,
creating a C++ class instance corresponding to the model using the shared object
code

We also want to be able to specify all the inputs and outputs including
adaptation, so that you can

* extract (and then potentially save) current adaptation parameters

* start a sampler with given adaptation parameters

Most of the time in our algorithm is spent computing gradients.
The faster that goes, the faster everything goes. If the adaptation
gets messed up, you'll have to intentionally slow everything down by
forcing a very small step size, which can deal with different situations.

The fact that interpreted R is snail slow isn't usually a problem if
you're just calling other functions that do the heavy lifting. But
you may find it prohibitive if you're in a tight loop doing a lot of
computation using R. But then if you have a good general technique
prototyped in R, it's usually pretty straightforward to move it to
something faster like C, C++, Fortran, or Java.

- Bob

[Krzysztof Sakrejda]

Thank you for taking the time to rewrite and explain where things are at.
This, combined with the optimization, is going to be amazing (even more
amazing than Stan already is). 

Krzysztof

[Michael Betancourt]

I think everyone should temper their expectations until

(a) we limit Stan models to DAGs or (b) we find some

magical algorithm that can compute marginal distributions

in closed form.

[Krzysztof Sakrejda]

On Wednesday, July 9, 2014 6:29:15 PM UTC-4, Michael Betancourt wrote:

I think everyone should temper their expectations until

(a) we limit Stan models to DAGs or (b) we find some

magical algorithm that can compute marginal distributions

in closed form.

Following Bob's notation, from Stan's point of view the distribution
p(alpha|beta,y) is going to be p(alpha|data)... Stan can compute that
density given data = {beta,y} as well as sample from it.  I don't think
we need p(alpha|y) if the goal is to have this be part of a Gibbs step. 

... what am I missing?

Krzysztof

[Bob Carpenter]

I don't see anything you're missing, so maybe Michael
will chime in.

I don't understand why Michael's talking about computing marginal
distributions in closed form --- we just have Bayes's rule and compute
everything up to a constant. I think that'll work for your
proposal.

If we did limit models to DAGs a la BUGS, then there's more
that we could infer from the graph structure, such as Markov
blankets, which would make doing a proper per-variable
Gibbs easy to implement. The fact that we don't do that is
why we don't have any discrete sampling.

And just to be clear about our goals, we're not planning on building these
Gibbs updates any time soon. And if you had to implement full Gibbs this way,
as in trying to code up BUGS in Stan, it'd be too much copying to be efficient
unless there were only a few blocks.

- Bob

[Sebastian Weber]

Hi!

Since you asked for features:

- a facility to simulate from a posterior would be nice, i.e. currently I loop over each draw of a posterior and fire up stan for each draw; the posterior should come either from a stanfit object or an array (I may want to mess around with the posterior before I simulate).

- some way to access the data from within the init function, i.e. like Ben pointed out; not so urgent, I am happy to live with the environment stuff I digged up in R.

- I saw you intend to parallelize things easily for users; here it would be great to keep things as flexible as possible (most people use some openmp, some MPI, and ...); concerning seeding I think it would be best to go with one stan seed which is the same for all chains, but different chain_ids; one more note: I keep reading that parallelization is not possible under windows which is not true, it is just different, i.e. you need to use "cl <- makeCluster(ncpu, type=ifelse(.Platform$OS.type == "windows", "PSOCK", "FORK"))"; but this PSOCK cluster needs to be handled with a bit of care (seeds, data exporting!); I use this with foreach which I find useful, but other will disagree here I guess

- consider to make the model chaching option configurable via the options interface from R so that users can set in their R profile if they like caching permanently on or not

- possibly consider to cache models in different compiled versions (fast/debug/NDEBUG) when they are requested in different flavors along the road

- I hope the current stan merge function stays, it is very useful as I run multiple chains on a cluster and each in a separate process which I merge at the end of the day

- Speeding up the monitor function would be great - I know you consider moving that bit into C++ to make it fast. If you could also write a function such that I can stick coda samples into it, that would even be better (sometimes I would like to do comparisons of Stan runs with JAGS/WinBUGS)

@ Bob: R is full of thousands of very convenient functions and I agree your closure type approach is programmatically much nicer (but the flow of my program does not go with it and after all it is just R - so who cares anyway).

Best,
Sebastian

[Sebastian Weber]

Hi!

One more thing: If possible, speed up the communication between R and Stan, i.e. I just had to rewrite a Stan model to expose fewer variables just because it took too long otherwise. I suspect that the issue is that Stan outputs all parameters while Rstan lets me select a subset and this subsetting is done on the R side. It seems to me that reading in large data chunks which come from Stan somehow takes a while (feels like stan is writing a csv file and rstan reads it in, which means rstan has to figure out the layout of the arrays).

Anything which speeds this stan<->R communication up and save me from copying my Stan models (with minor modifications then) around would be great.

Best,
Sebastian

[Krzysztof Sakrejda]

On Wednesday, July 9, 2014 11:25:38 PM UTC-4, Bob Carpenter wrote:

I don't see anything you're missing, so maybe Michael
will chime in.

Ok, this is something I'm doing in an applied project so I would like to
hear if there are serious problems with it.

I don't understand why Michael's talking about computing marginal
distributions in closed form --- we just have Bayes's rule and compute
everything up to a constant.  I think that'll work for your
proposal.  

If we did limit models to DAGs a la BUGS, then there's more
that we could infer from the graph structure, such as Markov
blankets, which would make doing a proper per-variable
Gibbs easy to implement.  The fact that we don't do that is
why we don't have any discrete sampling.

And just to be clear about our goals, we're not planning on building these
Gibbs updates any time soon. And if you had to implement full Gibbs this way,
as in trying to code up BUGS in Stan, it'd be too much copying to be efficient
unless there were only a few blocks.

Yeah, a lot of copying.  I did a proof-of-concept a while ago which I'm pasting below.
Note that it currently does not work b/c I've changed the package code (and it requries
a library I have on github, which has also changed).  I am currently updating my examples
and I will link them here it in the next week or so to make it run-able.  For now it's an ugly mock-up...
but at least it makes my goals clear.  I think the ability to do a few big blocks will save me
and maybe others from writing NUTS by hand or writing the umpteenth version of slow code for doing an HLM.

library(Rux)
library(MASS)
library(rstan)


## From total marked population size (N), times of surveys (tos),
## survival probability (phi), and recapture probability (p),
## simulate times of ``births'', followed by the survival process
## (time of alive matrix, toaM), and the recapture process
## (time of recpature matrix, torM).  From these, calculate the
## list of times of recapture (tor) and the vector of times of
## death (tod).
tos <- 1:10
N <- 200
phi <- rep(.8,max(tos))
p <- 0.7
tbs <- rbinom(n=N, size=5, prob=.3)+1
toaM <- matrix(0,nrow=N, ncol=length(tos))
torM <- matrix(0,nrow=N, ncol=length(tos))

survive <- 0
die <- 0
recapture <- 0
miss <- 0

for ( i in 1:N ) {
  toaM[i,tbs[i]] <- 1
  torM[i,tbs[i]] <- 1
  for (j in (tbs[i]+1):max(tos) ) {
    toaM[i,j] <- toaM[i,j-1] *
      sample(x=c(0,1),size=1,replace=TRUE, prob=c(1-phi[j], phi[j]))
    if ((toaM[i,j] == 1) && (toaM[i,j-1] == 1)) survive <- survive + 1
    if ((toaM[i,j] == 0) && (toaM[i,j-1] == 1)) die <- die + 1

    torM[i,j] <- toaM[i,j] *
      sample(x=c(0,1),size=1,replace=TRUE, prob=c(1-p, p))
    if ((toaM[i,j] == 1) && (torM[i,j] == 1)) recapture <- recapture + 1
    if ((toaM[i,j] == 1) && (torM[i,j] == 0)) miss <- miss + 1

  }

}
tor <- apply(X=torM, MARGIN=1, FUN=function(x) {which(x==1)})
tod <- apply(X=toaM, MARGIN=1, FUN=function(x) {max(which(x==1))+1})
                                               

## Calculate realized survival and recapture probabilities.
true_phi <- survive/(survive + die)
true_p <-   recapture/(recapture + miss)


## Create matrices of (bad) initial values for the parameters of the
## survival (PHI) and recapture (P) processes.  These matrices are
## N x max(tos)*5 by default.
P <- matrix(data=0.6, nrow=N, ncol=max(tos)*5)
P[,(max(tos)+1):(max(tos)*5)] <- 0
PHI = 0.5

## Create a handle to the likelihood calculation and discrete slice
## sampler for time of death state.
o <- new("recapture_handle",
  times_of_surveys = tos,
  times_of_recaptures = tor,
  times_of_deaths = tod,
  known_deaths = rep(FALSE, N),
  PHI = PHI,
  P = P
)

## Calculate the initial likelihood, the time taken for this calculation
## should scale with N * K * 5.  
ll_calculation_time <- system.time(expr=o$log_likelihood)


## Create vectors to hold output from MCMC.
N_it <- 1000
o_accept <- vector(mode="numeric", length=N_it)
o_PHI <- vector(mode="numeric", length=N_it)
o_P <- vector(mode="numeric", length=N_it)
o_lpd <- vector(mode="numeric", length=N_it)

## Stan model for simulating from PHI/P conditional on known
## birth, death, and recapture times for n tags.
stan_code <- '
  data {
    int<lower=0> n_tags;
    int<lower=0> n_occasions;
    int<lower=0> tb[n_tags];
    int<lower=0> td[n_tags];
    int<lower=0, upper=1> recaptures[n_tags, n_occasions];

  }

  // Declare parameters for the linear models.
  parameters {
    real beta;
    real gamma;
  }

  // Rubber meets the road for linear model specification.
  transformed parameters {
    real PHI_logit[n_tags, 5*n_occasions];
    real P_logit[n_tags, 5*n_occasions];
    for ( tag in 1:n_tags ) { // indexing covers PHI/P
      for ( t in tb[tag]:(td[tag]-1) ) {
        PHI_logit[tag,t] <- beta;
        P_logit[tag,t] <- gamma;
      }
    }
  }
 model {

    // Priors for linear model coefficients.
    beta ~ cauchy(0,2.5);
    gamma ~ cauchy(0,2.5);

    // Mark-recapture likelihood.
    for ( tag in 1:n_tags ) {
      for ( t in tb[tag]:(td[tag]-2) ) { // SURVIVES
        1 ~ bernoulli_logit(PHI_logit[tag,t]);
      }
      0 ~ bernoulli_logit(PHI_logit[tag,td[tag]-1]);  // DIES

      for ( t in (tb[tag]+1):min(n_occasions,td[tag]-1) ) {  // AVAILABLE
        recaptures[tag,t] ~ bernoulli_logit(P_logit[tag,t]);
      }
    }
  }

'

n_tags <- dim(o$recaptures)[1]
n_occ <- dim(o$recaptures)[2]
stan_fit <- NULL
stan_iter <- 200

stan_timing <- system.time({
for ( it in 1:N_it ) {
  print(it)
  ## Sample state:
  tds <- c(o$set_deaths(o$td))

  ## Sample PHI/P matrices:
  stan_fit <- stan(
    model_code=stan_code,
    data=list(
      n_tags = n_tags,
      n_occasions = n_occ,
      tb = tbs,
      td = tds,
      recaptures = o$recaptures
    ), fit=stan_fit,
    iter=stan_iter, chains=1, warmup=stan_iter/2,
    pars=c('beta','gamma')
  )

  ## Make matrices and re-assign them:
  stan_samples <- extract(stan_fit, permuted=TRUE)
  PHImat <- matrix(
    data=plogis(stan_samples$beta[stan_iter/2]),
    nrow=o$N, ncol=o$K*5, byrow=TRUE)
  Pmat <- matrix(
    data=plogis(stan_samples$gamma[stan_iter/2]),
    nrow=o$N, ncol=o$K*5, byrow=TRUE)
  Pmat[,(o$K+1):ncol(Pmat)] <- 0
  o$PHI <- PHImat
  o$P   <- Pmat

  ## Record for posterity
  o_PHI[it] = unique(c(o$PHI))
  o_P[it]   = max(unique(c(o$P)))
  o_lpd[it] = o$log_likelihood
}
})

plot(o_PHI, pch='*')
plot(o_lpd, pch='*')
plot(o_PHI, o_lpd, pch = '*')
plot(o_PHI, o_P, pch = '*')
abline(h=recapture/(recapture+miss))
abline(v=survive/(survive+die))

 

[Andrew Gelman]

We’re never gonna limit Stan models to Dags. I think what we’re going to build sometimes is a wrapper for Stan where the user can supply a Dag and also say which nodes are observed, and then the wrapper will write the appropriate Stan model.

On Jul 10, 2014, at 5:25 AM, Bob Carpenter <ca...@alias-i.com> wrote:
>
> If we did limit models to DAGs a la BUGS, then there's more
> that we could infer from the graph structure, such as Markov
> blankets, which would make doing a proper per-variable
> Gibbs easy to implement. The fact that we don't do that is
> why we don't have any discrete sampling.
>
>

[Bob Carpenter]

Right --- we're not going to limit them everywhere.

I just meant that for the special case of directed acyclic
graphical models, we could exploit the graph structure.

The real trick in Andrew's proposal would be handling
missing data a la BUGS --- we'd need to break the parameters
apart into observed (data) and "missing" (parameters), but then
I don't know how the interfaces on the outside could work
without major overhauls.

- Bob

[Andrew Gelman]

I hadn’t been thinking about that particular difficulty. It could make sense to do it using our current hack-structure for missing data. But I agree that if we could handle missing data in a natural way, that would be cool.

[Bob Carpenter]

Nothing is being written out to .csv format and
read back in. That would be super slow.
Everything communicates directly in memory through Rcpp.

Restricting the number of parameters output more efficiently
from the C++ Stan side is on our to-do list. But R should
already be able to do it pretty efficiently (in the sense
that the time should be negligible compared to computing
gradients each iteration).

I don't know where the hotspot in RStan is for this, but I'd
guess it's in reshaping everything to go into the fit
object. One approach we could take would be to do more work
lazily --- I've been thinking about that for the new interfaces.

- Bob

[Michael Betancourt]

Can you guys clarify something for me?  Are we trying to expose the

full conditional density,

p(alpha | beta, x)

or just something proportional to it,

~ p(alpha | beta, x) ~ p(alpha, beta) | beta=const?

The former is what’s hard and the cause of my concern.  If it’s just

the latter that’s necessary then it’s no problem.  Apologies if I was

mistaken.

[Jiqiang Guo]

RStan saves the samples directly to an R vector.  So I don't think

there is much to be improved.  Also there is not reshaping for the final fit object, so there is no extra overload. 

But if rstan is unexpected slow for a model with not many parameters, we should take a look and see if there are places that can be proved.  So please let us know such specific cases. 

Jiqiang 

[Krzysztof Sakrejda]

On Thursday, July 10, 2014 5:46:23 PM UTC-4, Michael Betancourt wrote:

Can you guys clarify something for me?  Are we trying to expose the

full conditional density,

p(alpha | beta, x)

or just something proportional to it,

~ p(alpha | beta, x) ~ p(alpha, beta) | beta=const?

The former is what’s hard and the cause of my concern.  

I think in all cases we can work with something proportional to p(alpha |beta, x).

The two steps are: a) simulate alpha from something proportional to p(alpha|beta,x),
let's say Stan can do this; and b) simulate beta from something proportional to 
p(beta|alpha,x), let's say something outside of Stan (say R, Python) can do this.
I don't think we ever need the full conditional, just something proportional.

I don't think the density proportional to p(alpha|beta,x) ever needs to be exposed from
Stan for this scheme to work. I just need to be able to read/write parameter/data efficiently.
It might be nice to get the density proportional to p(alpha|beta,x) but it is not necessary.

 

If it’s just

the latter that’s necessary then it’s no problem.  Apologies if I was

mistaken.

Thank you for chiming in, I've learned a lot from the Stan team's development
list and responses user questions and I appreciate the time you take to
resolve these questions.

Krzysztof
 

[Sebastian Weber]

It would be more than cool! Handling missing data automatically is still one of the real nice things BUGS has over Stan, i.e. not forcing users to do the index juggling would be nice. I guess this would also enable things like easy simulation from a model and automatic WAIC calculations - even better!

The approach of generating from a DAG description file a Stan model sounds perfect (because I can then alter the file generated and hence get the benefits of DAGs and are still not limited by them).

Is that for Stan 4.0?

Best,
Sebastian

[Michael Betancourt]

The two steps are: a) simulate alpha from something proportional to p(alpha|beta,x), 
let's say Stan can do this; and b) simulate beta from something proportional to  
p(beta|alpha,x), let's say something outside of Stan (say R, Python) can do this.
I don't think we ever need the full conditional, just something proportional.

I don't think the density proportional to p(alpha|beta,x) ever needs to be exposed from
Stan for this scheme to work. I just need to be able to read/write parameter/data efficiently.
It might be nice to get the density proportional to p(alpha|beta,x) but it is not necessary.

Now I’m hesitant again.  In order to, say, sample a conditional from within Stan we’d have

to compute the joint (proportional to the conditional).  No external program will be able to

do the same without accessing the density, at the very least.

[Sebastian Weber]

Hi!

Maybe I overlooked it, but a function in rstan like this one:

check_sampler <- function(fit, ...) {
    sp <- get_sampler_params(fit, inc_warmup=FALSE)
    spa <- do.call(abind, c(list(along=-1), sp))  
    spa <- aperm(spa, c(2,1,3))
    monitor(spa, warmup=0, ...)
}

would be nice in order to check that the sample did fine. My version here needs the abind package.

Best,
Sebastian

--

[Krzysztof Sakrejda]

I think we're still on the same page: I expect the outside program to also implement the
density so it can carry out the sampling.  Maybe we're talking past each other so when
I get a more transparent working example I'll post it here.

Krzysztof

[Ben Goodrich]

On Friday, July 11, 2014 5:00:17 AM UTC-4, Sebastian Weber wrote:

Maybe I overlooked it, but a function in rstan like this one:

check_sampler <- function(fit, ...) {
    sp <- get_sampler_params(fit, inc_warmup=FALSE)
    spa <- do.call(abind, c(list(along=-1), sp))  
    spa <- aperm(spa, c(2,1,3))
    monitor(spa, warmup=0, ...)
}

would be nice in order to check that the sample did fine. My version here needs the abind package.

You can check whether fit@mode is zero to verify that there are samples. Otherwise, I think

monitor(as.array(fit), warmup = 0)

would do pretty much the same thing.

Ben

[Sebastian Weber]

Hi Ben!

If I run as.array(fit), then I the sampler parameters are not part of the output (also the permuted argument seems broken):

R3.0.2> aa <- as.array(f$fit, permuted=FALSE)
Error in .local(object, ...) :
  formal argument "permuted" matched by multiple actual arguments

Moreover, I noticed that monitor does not like everything I give it:
eight <- stan_demo("eight_schools", chains=1)
print(eight)
check_sampler(eight)

The array is correctly formatted, but monitor in check_sampler gets confused for reasons I don't get.

Best,

Sebastian

[Andrew Gelman]

S

I don’t know when this will be done—they key difficulty is that someone would have to do it!—but I don’t think it has to wait on any other aspects of Stan.  At least as I’ve thought about it, we’re talking about a wrapper for Stan rather than a change in Stan itself.

A

[Ben Goodrich]

On Friday, July 11, 2014 10:33:28 AM UTC-4, Sebastian Weber wrote:

Hi Ben!

If I run as.array(fit), then I the sampler parameters are not part of the output (also the permuted argument seems broken):

R3.0.2> aa <- as.array(f$fit, permuted=FALSE)
Error in .local(object, ...) :
  formal argument "permuted" matched by multiple actual arguments

The as.array() method hard-codes permuted=FALSE, so it generates an error if you try to pass it

> rstan:::as.array.stanfit
function (x, ...)
{
    if (x@mode != 0)
        return(numeric(0))
    out <- extract(x, permuted = FALSE, inc_warmup = FALSE, ...)
    return(out)
}
<environment: namespace:rstan>

but I see now that doesn't produce the stuff that you want.

Moreover, I noticed that monitor does not like everything I give it:

eight <- stan_demo("eight_schools", chains=1)
print(eight)
check_sampler(eight)

The array is correctly formatted, but monitor in check_sampler gets confused for reasons I don't get.

In many cases, n_divergent is always zero and thus it throws an error when trying to inflate the variance. We'll fix that, but for now you could do

check_sampler <- function(fit, ...) {
  sp <- get_sampler_params(fit, inc_warmup=FALSE)
  spa <- do.call(abind, c(list(along=-1), sp)) 
  spa <- aperm(spa, c(2,1,3))

  spa <- spa[,,apply(spa, 3, sd) > 0, drop = FALSE]
  monitor(spa, warmup=0, ...)
}

Ben

[Bob Carpenter]

It's not on our short-term schedule and we don't have
much of a long-term schedule.

This could actually be a nice standalone project. It depends
on Stan, but nothing in Stan should need to be changed to
get it all going. It'd be a wrapper.

We'd need someone who could build a parser and then do the
code generation, but otherwise, this shouldn't be too difficult.

- Bob

[Bob Carpenter]

I missed Andrew's response --- what he said!

- Bob

[Richard McElreath]

I'm not sure what extra is required for DAG desires, but my R package "rethinking" (github.com/rmcelreath/rethinking) has a function map2stan() that takes formula input like this:

f <- alist(

    y ~ dnorm( mu , sigma ),

    mu <- a + b*x,

    a ~ dnorm(0,10),

    b ~ dnorm(0,1),

    sigma ~ dexp(1)

)

and it produces Stan code like:

data{

    int<lower=1> N;

    real y[N];

    real x[N];

}

parameters{

    real a;

    real b;

    real<lower=0> sigma;

}

model{

    vector[N] mu;

    sigma ~ exponential( 1 );

    b ~ normal( 0 , 1 );

    a ~ normal( 0 , 10 );

    for ( i in 1:N ) {

        mu[i] <- a + b * x[i];

    }

    y ~ normal( mu , sigma );

}

It wouldn't be hard to tag nodes in this interface. It support a few posterior prediction functions:

* link() to get values of linear model(s) at each sample

* sim() to get samples from outcome at each sample

* postcheck() just computes posterior predictions for each case in the sample and plots them

* WAIC() to calculate WAIC

Anyway, if I can get a better idea what the desired functionality is, it might be easier to add it to my package than to write another front end. Or you folks can just fork my work, although I imagine learning my parsing algorithm (it uses templates to map distributions from R to Stan) would be as much work as writing your own.

[Bob Carpenter]

Cool. We totally have to check that out.

For functionality:

1. Basically, the whole BUGS/JAGS language. Ideally,
something that would translate models in both of those
languages!

2. Something that would take in the data for the model,
figure out which values are "missing" and create a model
where the missing data are considered parameters and then
glue everything back together. I'm not even sure how
BUGS does this in terms of returns, since the missing data
can be part of bigger structures (at least in some cases).

3. Something that could compute the Markov blanket for
discrete parameters so that we could do discrete sampling
with Gibbs. This is a wide-open research area in terms of
how this will interact with HMC and how adaptation would work.

4. Whatever Andrew wants for model comparison, data simulation
and posterior predictive checks. Looks like you're a long way
down that road already.


For compatibility:

1. We want everything we build for Stan to work in
CmdStan and PyStan as well as RStan. So we couldn't
have dependencies to R or anything else that's not
BSD-compatibly licensed (specifically, no GPL).

- Bob

P.S. We're also interested in something similar for the
R expression language --- I know you've worked on that before,
too. We're working on implementing uncertainty estimation for
MLE and also doing marginal maximum likelihood. There are still
some issues in how we take uncertainty back from the unconstrained
parameters to the constrained --- the best idea so far is to
use the normal approximation to generate samples around the mode
then transform back to the constrained space.