stan-dev / stan-dev/rstan

read_stan_csv parameter naming scheme differs from sampling output

Open
#733 2 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

Dominant language
R
Stars
1.1k
Forks
266
Avg merge
2h 56m
Merged PRs (30d)
1

Description

Summary:

The naming scheme used to reconstruct the stanfit object in read_stan_csv is different than the original rstan::sampling output

shredder::rat_data

$N
[1] 30

$T
[1] 5

$x
[1]  8 15 22 29 36

$y
# A tibble: 30 x 5
    day8 day15 day22 day29 day36
   <int> <int> <int> <int> <int>
 1   151   199   246   283   320
 2   145   199   249   293   354
 3   147   214   263   312   328
 4   155   200   237   272   297
 5   135   188   230   280   323
 6   159   210   252   298   331
 7   141   189   231   275   305
 8   159   201   248   297   338
 9   177   236   285   350   376
10   134   182   220   260   296
11   160   208   261   313   352
12   143   188   220   273   314
13   154   200   244   289   325
14   171   221   270   326   358
15   163   216   242   281   312
16   160   207   248   288   324
17   142   187   234   280   316
18   156   203   243   283   317
19   157   212   259   307   336
20   152   203   246   286   321
21   154   205   253   298   334
22   139   190   225   267   302
23   146   191   229   272   302
24   157   211   250   285   323
25   132   185   237   286   331
26   160   207   257   303   345
27   169   216   261   295   333
28   157   205   248   289   316
29   137   180   219   258   291
30   153   200   244   286   324

$xbar
[1] 22



system.file('rats.stan',package='shredder')

data {
  int<lower=0> N;
  int<lower=0> T;
  real x[T];
  real y[N,T];
  real xbar;
}
parameters {
  real alpha[N];
  real beta[N];

  real mu_alpha;
  real mu_beta;          // beta.c in original bugs model

  real<lower=0> sigmasq_y;
  real<lower=0> sigmasq_alpha;
  real<lower=0> sigmasq_beta;
}
transformed parameters {
  real<lower=0> sigma_y;       // sigma in original bugs model
  real<lower=0> sigma_alpha;
  real<lower=0> sigma_beta;

  sigma_y = sqrt(sigmasq_y);
  sigma_alpha = sqrt(sigmasq_alpha);
  sigma_beta = sqrt(sigmasq_beta);
}
model {
  mu_alpha ~ normal(0, 100);
  mu_beta ~ normal(0, 100);
  sigmasq_y ~ inv_gamma(0.001, 0.001);
  sigmasq_alpha ~ inv_gamma(0.001, 0.001);
  sigmasq_beta ~ inv_gamma(0.001, 0.001);
  alpha ~ normal(mu_alpha, sigma_alpha); // vectorized
  beta ~ normal(mu_beta, sigma_beta);  // vectorized
  for (n in 1:N)
    for (t in 1:T) 
      y[n,t] ~ normal(alpha[n] + beta[n] * (x[t] - xbar), sigma_y);

}
generated quantities {
  real alpha0;
  alpha0 = mu_alpha - xbar * mu_beta;
}


library(shredder)
library(rstan)
#> Loading required package: StanHeaders
#> Loading required package: ggplot2
#> rstan (Version 2.19.2, 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)

junk <- capture.output({ #verbose doesnt work in reprex
  fit <- rstan::stan(
    file  = system.file('rats.stan',package='shredder'),
    data  = shredder::rat_data,verbose = FALSE
  )  
})
#> Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
#> http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: The largest R-hat is 1.05, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#r-hat
#> 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
#> http://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
#> http://mc-stan.org/misc/warnings.html#tail-ess

m <- stan_model(system.file('rats.stan',package='shredder'))

td <- file.path(tempdir(),'test')
dir.create(td)

junk <- capture.output({ #verbose doesnt work in reprex
f <- sampling(m, data = shredder::rat_data,sample_file = file.path(td,'samp'),verbose = FALSE)
})

from_csv <- rstan::read_stan_csv(list.files(td,full.names = TRUE))

class(from_csv)
#> [1] "stanfit"
#> attr(,"package")
#> [1] "rstan"

class(fit)
#> [1] "stanfit"
#> attr(,"package")
#> [1] "rstan"

names(from_csv@sim$samples[[1]])
#>  [1] "alpha.1"       "alpha.2"       "alpha.3"       "alpha.4"      
#>  [5] "alpha.5"       "alpha.6"       "alpha.7"       "alpha.8"      
#>  [9] "alpha.9"       "alpha.10"      "alpha.11"      "alpha.12"     
#> [13] "alpha.13"      "alpha.14"      "alpha.15"      "alpha.16"     
#> [17] "alpha.17"      "alpha.18"      "alpha.19"      "alpha.20"     
#> [21] "alpha.21"      "alpha.22"      "alpha.23"      "alpha.24"     
#> [25] "alpha.25"      "alpha.26"      "alpha.27"      "alpha.28"     
#> [29] "alpha.29"      "alpha.30"      "beta.1"        "beta.2"       
#> [33] "beta.3"        "beta.4"        "beta.5"        "beta.6"       
#> [37] "beta.7"        "beta.8"        "beta.9"        "beta.10"      
#> [41] "beta.11"       "beta.12"       "beta.13"       "beta.14"      
#> [45] "beta.15"       "beta.16"       "beta.17"       "beta.18"      
#> [49] "beta.19"       "beta.20"       "beta.21"       "beta.22"      
#> [53] "beta.23"       "beta.24"       "beta.25"       "beta.26"      
#> [57] "beta.27"       "beta.28"       "beta.29"       "beta.30"      
#> [61] "mu_alpha"      "mu_beta"       "sigmasq_y"     "sigmasq_alpha"
#> [65] "sigmasq_beta"  "sigma_y"       "sigma_alpha"   "sigma_beta"   
#> [69] "alpha0"        "lp__"

names(fit@sim$samples[[1]])
#>  [1] "alpha[1]"      "alpha[2]"      "alpha[3]"      "alpha[4]"     
#>  [5] "alpha[5]"      "alpha[6]"      "alpha[7]"      "alpha[8]"     
#>  [9] "alpha[9]"      "alpha[10]"     "alpha[11]"     "alpha[12]"    
#> [13] "alpha[13]"     "alpha[14]"     "alpha[15]"     "alpha[16]"    
#> [17] "alpha[17]"     "alpha[18]"     "alpha[19]"     "alpha[20]"    
#> [21] "alpha[21]"     "alpha[22]"     "alpha[23]"     "alpha[24]"    
#> [25] "alpha[25]"     "alpha[26]"     "alpha[27]"     "alpha[28]"    
#> [29] "alpha[29]"     "alpha[30]"     "beta[1]"       "beta[2]"      
#> [33] "beta[3]"       "beta[4]"       "beta[5]"       "beta[6]"      
#> [37] "beta[7]"       "beta[8]"       "beta[9]"       "beta[10]"     
#> [41] "beta[11]"      "beta[12]"      "beta[13]"      "beta[14]"     
#> [45] "beta[15]"      "beta[16]"      "beta[17]"      "beta[18]"     
#> [49] "beta[19]"      "beta[20]"      "beta[21]"      "beta[22]"     
#> [53] "beta[23]"      "beta[24]"      "beta[25]"      "beta[26]"     
#> [57] "beta[27]"      "beta[28]"      "beta[29]"      "beta[30]"     
#> [61] "mu_alpha"      "mu_beta"       "sigmasq_y"     "sigmasq_alpha"
#> [65] "sigmasq_beta"  "sigma_y"       "sigma_alpha"   "sigma_beta"   
#> [69] "alpha0"        "lp__"

R.version.string
#> [1] "R version 3.6.1 (2019-07-05)"

packageVersion("rstan")
#> [1] '2.19.2'

details::details(sessioninfo::session_info(),summary = 'R Session Info')
R Session Info

─ Session info ──────────────────────────────────────────────────────────
 setting  value                       
 version  R version 3.6.1 (2019-07-05)
 os       macOS Mojave 10.14.5        
 system   x86_64, darwin15.6.0        
 ui       X11                         
 language (EN)                        
 collate  en_US.UTF-8                 
 ctype    en_US.UTF-8                 
 tz       America/New_York            
 date     2020-01-17                  

─ Packages ──────────────────────────────────────────────────────────────
 package     * version    date       lib source                         
 assertthat    0.2.1      2019-03-21 [1] CRAN (R 3.6.0)                 
 backports     1.1.5      2019-10-02 [1] CRAN (R 3.6.0)                 
 callr         3.3.2      2019-09-22 [1] CRAN (R 3.6.0)                 
 checkmate     1.9.4      2019-07-04 [1] CRAN (R 3.6.0)                 
 cli           2.0.1      2020-01-08 [1] CRAN (R 3.6.0)                 
 clipr         0.7.0      2019-07-23 [1] CRAN (R 3.6.0)                 
 codetools     0.2-16     2018-12-24 [1] CRAN (R 3.6.1)                 
 colorspace    1.4-1      2019-03-18 [1] CRAN (R 3.6.0)                 
 crayon        1.3.4      2017-09-16 [1] CRAN (R 3.6.0)                 
 desc          1.2.0      2019-12-01 [1] Github (r-lib/desc@61205f6)    
 details       0.2.1      2020-01-12 [1] local                          
 digest        0.6.23     2019-11-23 [1] CRAN (R 3.6.0)                 
 dplyr         0.8.3      2019-07-04 [1] CRAN (R 3.6.0)                 
 evaluate      0.14       2019-05-28 [1] CRAN (R 3.6.0)                 
 fansi         0.4.1      2020-01-08 [1] CRAN (R 3.6.0)                 
 ggplot2     * 3.2.1      2019-08-10 [1] CRAN (R 3.6.0)                 
 glue          1.3.1.9000 2020-01-07 [1] Github (tidyverse/glue@b9ffe6c)
 gridExtra     2.3        2017-09-09 [1] CRAN (R 3.6.0)                 
 gtable        0.3.0      2019-03-25 [1] CRAN (R 3.6.0)                 
 highr         0.8        2019-03-20 [1] CRAN (R 3.6.0)                 
 htmltools     0.4.0      2019-10-04 [1] CRAN (R 3.6.0)                 
 httr          1.4.1      2019-08-05 [1] CRAN (R 3.6.0)                 
 inline        0.3.15     2018-05-18 [1] CRAN (R 3.6.0)                 
 knitr         1.25       2019-09-18 [1] CRAN (R 3.6.0)                 
 lazyeval      0.2.2      2019-03-15 [1] CRAN (R 3.6.0)                 
 lifecycle     0.1.0      2019-08-01 [1] CRAN (R 3.6.0)                 
 loo           2.1.0      2019-03-13 [1] CRAN (R 3.6.0)                 
 magrittr      1.5        2014-11-22 [1] CRAN (R 3.6.0)                 
 matrixStats   0.54.0     2018-07-23 [1] CRAN (R 3.6.0)                 
 munsell       0.5.0      2018-06-12 [1] CRAN (R 3.6.0)                 
 pillar        1.4.3      2019-12-20 [1] CRAN (R 3.6.0)                 
 pkgbuild      1.0.6      2019-10-09 [1] CRAN (R 3.6.1)                 
 pkgconfig     2.0.3      2019-09-22 [1] CRAN (R 3.6.0)                 
 png           0.1-7      2013-12-03 [1] CRAN (R 3.6.0)                 
 prettyunits   1.0.2      2015-07-13 [1] CRAN (R 3.6.0)                 
 processx      3.4.1      2019-07-18 [1] CRAN (R 3.6.0)                 
 ps            1.3.0      2018-12-21 [1] CRAN (R 3.6.0)                 
 purrr         0.3.3      2019-10-18 [1] CRAN (R 3.6.0)                 
 R6            2.4.1      2019-11-12 [1] CRAN (R 3.6.0)                 
 Rcpp          1.0.3      2019-11-08 [1] CRAN (R 3.6.1)                 
 rematch2      2.1.0      2019-07-11 [1] CRAN (R 3.6.0)                 
 rlang         0.4.2      2019-11-23 [1] CRAN (R 3.6.0)                 
 rmarkdown     2.0        2019-12-12 [1] CRAN (R 3.6.0)                 
 rprojroot     1.3-2      2018-01-03 [1] CRAN (R 3.6.0)                 
 rstan       * 2.19.2     2019-07-09 [1] CRAN (R 3.6.0)                 
 scales        1.1.0      2019-11-18 [1] CRAN (R 3.6.0)                 
 sessioninfo   1.1.1      2018-11-05 [1] CRAN (R 3.6.0)                 
 shredder    * 0.1.1      2019-12-20 [1] local                          
 StanHeaders * 2.18.1-10  2019-06-14 [1] CRAN (R 3.6.0)                 
 stringi       1.4.3      2019-03-12 [1] CRAN (R 3.6.0)                 
 stringr       1.4.0      2019-02-10 [1] CRAN (R 3.6.0)                 
 tibble        2.1.3      2019-06-06 [1] CRAN (R 3.6.0)                 
 tidyselect    0.2.5      2018-10-11 [1] CRAN (R 3.6.0)                 
 withr         2.1.2      2018-03-15 [1] CRAN (R 3.6.0)                 
 xfun          0.10       2019-10-01 [1] CRAN (R 3.6.0)                 
 xml2          1.2.2      2019-08-09 [1] CRAN (R 3.6.0)                 
 yaml          2.2.0      2018-07-25 [1] CRAN (R 3.6.0)                 

[1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library

Created on 2020-01-17 by the reprex package (v0.3.0)

Contributor guide

No contributing guide indexed for this repository

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

Research direction

Start with rstan::read_stan_csv and sampling, reproducing the discrepancy with the provided rats.stan example and shredder::rat_data. Compare the parameter names in the reconstructed and sampled stanfit objects; done means read_stan_csv uses the same naming scheme as sampling.

Written by the indexing model from the issue text.

Assessment

Tech stack
r
Domain
tooling
Issue type
Bug
Difficulty
3/5
Estimated time
1-2 days
Activity status
Stale
Clarity
Mostly clear
Newbie friendliness
42/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.