Various observations on rstan, cmdstanr, pystan and cmdstanpy after teaching with all in parallel

[amas0]

@robertgrant

In general, I am annoyed by Python packages that are called x for pip install and then when you go to import them, they are suddenly called y. Like pystan. Maybe that’s just me, but every time I explain to learners that they have to look out for things like that, it deflates their enthusiasm a little.

@ahartikainen

… I feel your pain, still I like the stan name :D

Just a thought on this – It seems possible to effectively alias pystan and stan in one installation in Python. The idea being you create a new package directory in the repo:

stan/
     __init__.py
     common.py
     fit.py
     model.py
     plugins.py
pystan/
     __init__.py
pyproject.toml
setup.cfg

And within the pystan/__init__.py file, you have:

from stan import *

Then everything in the stan namespace is copied over to the pystan namespace as well. Only additional thing you need to do is tell the build instructions to install both the stan and pystan packages, which is accomplished in the pyproject.toml:

packages = [
    { include = "stan" },
    { include = "pystan" },
]

I tested this out locally and it seems to work without issue. Could alleviate some of the issues that seem to come up when people install pystan and can’t import it with the same name. (Additionally, it looks like @ariddell has the stan package on pypi, so in theory you could have the build instructions push the package to both the pystan and stan pypi repositories.

I know the idea of having two ways to install/import the same package is kind of insane, but it’s a thought.

[amas0]

To continue hijacking this thread.

@robertgrant

It’s a shame about Jupyter notebooks, pystan and nest-asyncio. I don’t know enough about the issue to comment. Maybe there’s a way to get around it.

@ahartikainen

This is kind of a problem in Jupyter, which means it should be fixed there.
One option would always call nest-asyncio when using pystan, but I see many potential bugs in this solution.

One option would be to have pystan do some introspection on import to check whether or not it is being run within a Jupyter notebook. For example, putting the following code snippet in stan/__init__.py would run nest-asyncio whenever the package is imported in a notebook:

try:
    ipython_shell = get_ipython().__class__.__name__
    if ipython_shell == 'ZMQInteractiveShell':
        import nest_asyncio
        nest_asyncio.apply()
except NameError:
    pass

I tried it out locally and it worked fine – in practice, if people already need to invoke nest-asyncio when running in Jupyter, this could be a fine band-aid to use. It really doesn’t look like a proper fix is forthcoming any time soon.

[robertgrant]

Thanks for the ideas.

@ahartikainen

Currently both pystan and cmdstanpy has similar functions in the background (but they handle numpy arrays too), so there could a user facing function to use them.

That would be really helpful for anyone who isn’t highly experienced with Python. I give students the data in CSV, .rds for R and .json for Python, but the real world isn’t so kind! And I know most researchers don’t want to think about intermediate data storage formats, they just want to hit the button and get Bayesian!

@ahartikainen

pystan (httpstan actually) caches the model so there is not really recompilation but reading the model from cache. So it is ok to call the build command in a for loop multiple times.

That’s great – I didn’t know that. I still prefer compiling a sampler that contains a likelihood and a prior, and after compilation will take data, inits, seed, etc as input, but I can live with it.

[Funko_Unko]

@robertgrant

I still prefer compiling a sampler that contains a likelihood and a prior, and after compilation will take data, inits, seed, etc as input, but I can live with it.

I also find the current pystan behavior absolutely confusing. I want to build the model, and then maybe set the data that induce the posterior, but I don’t see how you can build a posterior.

[defjaf]

@Funko_Unko

I also find the current pystan behavior absolutely confusing. I want to build the model, and then maybe set the data that induce the posterior, but I don’t see how you can build a posterior.

I completely agree.

@ahartikainen

This was selected after long conversation (model + data → posterior → sample posterior)

Could you point to this long conversation? I would like to be convinced that this makes any sense at all. (I have commented elsewhere that it is a pain point for upgrading the jupyter/stan %magic functions which operate on a code block for compilation, and therefore don’t naturally have access to the data.)

[ahartikainen]

RStan3 and PyStan3 Interface

We had an unusually productive discussion today on Stan 3, mostly from the perspective of RStan. You can install a somewhat working version of the provisional rstan3 package from GitHub with devtools::install_github("stan-dev/rstan", ref = "develop", subdir = "rstan3") which sits on top of rstan 2.15.1 from CRAN. Much later I will make a branch that can be merged into rstan and then rstan3 will be deleted. I also updated the wiki but if you want to follow what rstan3 is doing, it is probably …

[defjaf]

Ok, well, thanks.

I guess I disagree with the outcome — users do, in fact, understand that the stan model program is separate from the specific data-set, and forcing those two things together as “build” breaks that mental model. The stan program is 100% independent of the values (but not the form) of the data, which is why stanc doesn’t need it.

In practice, with caching, this isn’t that big a deal in many use cases (but not all: the aforementioned Jupyter %stan magic just does make more sense with separate compilation and instantiation with data).

Moreover, I think it is not always a good pattern for a wrapper to ignore the implementation details, since at the very least it’s a good idea for users to understand what is going on in case things go wrong…

But this ship appears to have sailed.

[mike-lawrence]

Just chiming in to point out that the {model + data} framework also kinda misfits in the context of SBC, which in turn is becoming more common/recommended for workflows.

[mitzimorris]

@robertgrant

rstan and cmdstanr accept a list or a filename for data. Cmdstanpy and pystan want only a dictionary. Writing and reading from JSON in Python is a massive pain, at least in my experience, so it would be great if pystan and cmdstanpy took data files as an argument, and also offered their own data-writing function that could handle numpy arrays and the various float lengths that might entail.

CmdStanPy most certainly accepts filenames for data, and does provide its own data-writing function - it’s in both the doc and the very first “hello world” example: “Hello, World” — CmdStanPy 0.9.76 documentation

The data can be specified either as a filepath or a Python dict; in this example, we use the example datafile bernoulli.data.json:

here’s the doc for sample inference method, and methods for optimization, variational inference, and standalone generated quantities do the same:

data (Optional [ Union [ Dict, str ]] ) – Values for all data variables in the model,
specified either as a dictionary with entries matching the data variables,
or as the path of a data file in JSON or Rdump format.

@robertgrant

rstan allows a single function stan() to rule them all. I’m not saying I like rstan::stan(), just that it’s good to have options for everyone to pick from, especially beginners.

if the goal is to understand Stan and use existing models, then perhaps it would be better to only teach rstanarm and/or brms. If you’re teaching Stan in order to help people develop their own models, then they will most likely struggle with compilation, and getting a model to compile will be a nice milestone.

@robertgrant

I know readthedocs is convenient, and I’m not criticising anyone because my documentation is the worst, but really the hello world examples there are not sufficient for beginners, and the layout of the API reference takes some getting used to.

agreed, we need much more in the way of documentation. there is a neglected subfolder of jupyter notebooks. another way to go would be to translate existing case studies to CmdStanPy and I would be more than happy to help folks do this.

thanks for the wonderful and extensive feedback!

[mitzimorris]

@robertgrant

Suppose you want 1000 warmup draws and then 1000 retained after that. rstan requires warmup=1000, iter=2000, while all the others will have some form of num_warmup=1000,num_sample=1000. That’s a potential trap for the unwary.

+100 as the emoji goes

@robertgrant : would love to hear how the installation process went. we know it’s a huge pain point - any insights would be much appreciated

also, does choice of Python or R correspond to any difference in overall level of programming ability, specifically, ability to manage input and output data files, organize work, and write general-purpose functions? or are there lots of students who want to get into this data science thing without being able to work from the shell and navigate the filesystem among both R and Python users?

[robertgrant]

@mitzimorris

CmdStanPy most certainly accepts filenames for data, and does provide its own data-writing function - it’s in both the doc and the very first “hello world” example

Those examples seem to me to assume that bernoulli.data.json already exists.

I can’t comment on the internal data-file-writing functionality as it’s not in the API. It appears to call python json because it throws an error on numpy arrays in just the same way. It would be nice to have it available to all users. (I don’t know if “expose” is strictly the right term for that.) And if it can’t deal with numpy, then it would be handy to give some advice to users at least, and better still to provide a simple writing function that takes everything via strings. (I find myself pining for C++ here, and that does not happen often)

[robertgrant]

@mitzimorris

@robertgrant : would love to hear how the installation process went. we know it’s a huge pain point - any insights would be much appreciated

Absolutely fine in all interfaces from my experience, which is either Mac or Linux.

@mitzimorris

also, does choice of Python or R correspond to any difference in overall level of programming ability, specifically, ability to manage input and output data files, organize work, and write general-purpose functions? or are there lots of students who want to get into this data science thing without being able to work from the shell and navigate the filesystem among both R and Python users?

I don’t have enough data to say. It’s an important question though. What I can say with some confidence is that, across various training, coaching and consulting jobs, my prejudices about R vs Python users are constantly proven wrong.

“…you guys think I don’t give you straight answers. You have to talk to these statisticians. They will not give you a direct answer on anything. We don’t know. See, when did you start with the first two options, it’s either straight up and straight down, right? Or a total V, or maybe it’s up with a plateau and we’re somewhere on the plateau. They don’t know.” — 5 April 2020, Governor Andrew Cuomo COVID-19 Briefing April 5: Hospitalizations Dropping in NY - Rev, 20:35.

[mitzimorris]

@robertgrant

I can’t comment on the internal data-file-writing functionality as it’s not in the API. It appears to call python json because it throws an error on numpy arrays in just the same way. It would be nice to have it available to all users. (I don’t know if “expose” is strictly the right term for that.)

this is a good idea. filed issue need utility function "stan_json" · Issue #410 · stan-dev/cmdstanpy · GitHub

if you could please add an example of data that causes json errors, that would be great, as the existing utility function jsondump has logic that’s supposed to handle this.

stan-dev/cmdstanpy/blob/90d86dc8c88d20dadd8b96ea816cee1bf1f2d006/cmdstanpy/utils.py#L369-L379

[robertgrant]

Thanks for prioritising this @mitzimorris. I see that write_stan_json is now in the 0.9.77 release. I don’t have the exact setup that caused problems anymore; I’d have to spin up another server with it, so I will test it out for sure in the fullness of time.

[ariddell]

@ahartikainen

@robertgrant

I see some effort to converge on the names of arguments bearing fruit. It would be nice to have more of this. Of course, its possible to have two options available, one deprecated, like warmup= and num_warmup=

We probably would want to use the same keywords in all interfaces (now pystan and cmdstanpy differ a bit)

I thought we had a solution to this: we use cmdstan’s argument names. What’s left to do? Are there changes that need to be made in pystan?

[robertgrant]

rstan does its own thing, at least as of this summer.

[mitzimorris]

names for CmdStanX wrappers for the number of warmup and sampling iterations are iter_warmup and iter_sampling, respectively.

this was to make it clear that these are the number of iterations. because it’s possible to thin the reported draws, the number of draws returned is not necessarily the same as the number of iterations, sometimes folks got confused by what was meant by num_samples, iter_sampling is (perhaps) more clear.

the difference between rstan and cmdstan et al remains - in rstan iter is the total number of iterations and warmup is the number of that total which are warmup iterations, so that the number of sampling iterations is iter - warmup and if warmup isn’t specified, then the number of warmup and sampling iterations are both iter / 2.

[robertgrant]

This may sound pedantic, but I mention it just because you say the rationale for iter_* is that it’s the number of iterations… it’s not. It’s the number of accepted proposals. There are more iterations (not counting leapfrogging) that lead to proposals that are rejected. It’s like a while loop instead of a for loop.

[hhau]

@robertgrant

This may sound pedantic, but I mention it just because you say the rationale for iter_* is that it’s the number of iterations… it’s not. It’s the number of accepted proposals.

No, this is not true. See:

library(rstan)
#> Loading required package: StanHeaders
#> rstan (Version 2.26.0.9000, GitRev: 2e1f913d3ca3)
#> For execution on a local, multicore CPU with excess RAM we recommend calling
#> options(mc.cores = parallel::detectCores()).
#> To avoid recompilation of unchanged Stan programs, we recommend calling
#> rstan_options(auto_write = TRUE)
library(magrittr)
#> 
#> Attaching package: 'magrittr'
#> The following object is masked from 'package:rstan':
#> 
#>     extract
res <- stan(
  model_code = 'data {} parameters{real x;} model{x ~ std_normal();}',
  chains = 1,
  control = list(adapt_delta = 0.5)
)
#> SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 9e-06 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.007 seconds (Warm-up)
#> Chain 1:                0.006 seconds (Sampling)
#> Chain 1:                0.013 seconds (Total)
#> Chain 1:
#> Warning: There were 78 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
samples <- rstan::extract(res, pars = 'x', permuted = FALSE)

# if these were all 'accepted' then this would be a vector of length 100, 
# as each 'sample' would be an 'accepted' propsal (which would 
# necessarily differ)  from the original point, and thus
samples %>%
  magrittr::extract(1 : 1000, ,) %>% 
  round(digits = 8) %>% # floating point error will do funny things otherwise
  unique() %>% 
  length()
#> [1] 230

^{Created on 2021-10-07 by the reprex package (v2.0.1)}

adapt_delta is the target acceptance rate, so you can make unique(samples) differ arbitrarily from 1000 (well, between 1 and 1000) by decreasing its value.

[hhau]

I had to convince myself that this was also true for cmdstan:

library(cmdstanr)
#> This is cmdstanr version 0.4.0.9000
#> - Online documentation and vignettes at mc-stan.org/cmdstanr
#> - CmdStan path set to: /Users/amanderson/.cmdstanr/cmdstan-2.26.0
#> - Use set_cmdstan_path() to change the path
#> 
#> A newer version of CmdStan is available. See ?install_cmdstan() to install it.
#> To disable this check set option or environment variable CMDSTANR_NO_VER_CHECK=TRUE.
library(magrittr)

cmdstan_model <- cmdstan_model(
  write_stan_file('data {} parameters{real x;} model{x ~ std_normal();}')
)

run_model <- cmdstan_model$sample(
  chains = 1,
  iter_warmup = 200,
  iter_sampling = 1000,
  adapt_delta = 0.5
)
#> Running MCMC with 1 chain...
#> 
#> Chain 1 Iteration:    1 / 1200 [  0%]  (Warmup) 
#> Chain 1 Iteration:  100 / 1200 [  8%]  (Warmup) 
#> Chain 1 Iteration:  200 / 1200 [ 16%]  (Warmup) 
#> Chain 1 Iteration:  201 / 1200 [ 16%]  (Sampling) 
#> Chain 1 Iteration:  300 / 1200 [ 25%]  (Sampling) 
#> Chain 1 Iteration:  400 / 1200 [ 33%]  (Sampling) 
#> Chain 1 Iteration:  500 / 1200 [ 41%]  (Sampling) 
#> Chain 1 Iteration:  600 / 1200 [ 50%]  (Sampling) 
#> Chain 1 Iteration:  700 / 1200 [ 58%]  (Sampling) 
#> Chain 1 Iteration:  800 / 1200 [ 66%]  (Sampling) 
#> Chain 1 Iteration:  900 / 1200 [ 75%]  (Sampling) 
#> Chain 1 Iteration: 1000 / 1200 [ 83%]  (Sampling) 
#> Chain 1 Iteration: 1100 / 1200 [ 91%]  (Sampling) 
#> Chain 1 Iteration: 1200 / 1200 [100%]  (Sampling) 
#> Chain 1 finished in 0.0 seconds.
samples_cmdstan <- run_model$draws('x')
samples_cmdstan %>%
  magrittr::extract(1 : 1000, ,) %>% 
  round(digits = 8) %>% # floating point error will do funny things otherwise
  unique() %>% 
  length()
#> [1] 589

^{Created on 2021-10-07 by the reprex package (v2.0.1)}