# SBC for cmdstanr

**URL:** <https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931>\
**Category:** General\
**Tags:** simulation-based-calibration, cmdstanr\
**Created:** [October 27, 2020, 12:48pm UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931 "2020-10-27T12:48:17Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![bnicenboim](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bnicenboim/32/21197_2.png) [@bnicenboim](https://discourse.mc-stan.org/u/bnicenboim)\
**Post date:** [October 27, 2020, 12:48pm UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/1 "2020-10-27T12:48:17Z")

</div>

Not an inquiry, and I wasn’t sure where to post it.

I wrote a function to calculate SBC for a cmdstanr model in [https://gist.github.com/bnicenboim/5796795acf164f0a0f163401e5975bfd](https://gist.github.com/bnicenboim/5796795acf164f0a0f163401e5975bfd).  
For now parallelization is not working ([https://github.com/stan-dev/cmdstanr/issues/326](https://github.com/stan-dev/cmdstanr/issues/326)). But maybe this is useful for someone else.

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [October 27, 2020, 7:30pm UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/2 "2020-10-27T19:30:11Z")

</div>

Cool, thanks! We have an open issue for adding SBC functionality to CmdStanR

> <https://github.com/stan-dev/cmdstanr/issues/167>
>
> Is your feature request related to a problem? Please describe.
> rstan has some nascent support for simulation-based calibration
> Describe the solution you'd like
> It...

but haven’t gotten to it yet. Any interest in contributing this to the package?

---

<div class="post-metadata">

**Author:** ![bnicenboim](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bnicenboim/32/21197_2.png) [@bnicenboim](https://discourse.mc-stan.org/u/bnicenboim)\
**Post date:** [October 27, 2020, 7:54pm UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/3 "2020-10-27T19:54:41Z")

</div>

Sure, I could do it. I thought you wanted to keep cmdstanr lightweight, and I IMHO it would make more sense to have another package where tools are implemented for both interfaces. I have been using the rstan:::print.sbc and rstan:::plot.sbc to check the output. Should I just copy the code to cmdstanr, or is there a different preferred way to see the output?

Another thing is that I couldn’t parallelize it ([https://github.com/stan-dev/cmdstanr/issues/326](https://github.com/stan-dev/cmdstanr/issues/326)), I couldn’t find a straightforward way to get if the model finished without warnings/errors or not, and maybe less importantly I have no clue on how to save an unfinished calibration (as rstan does).

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [October 27, 2020, 7:56pm UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/4 "2020-10-27T19:56:45Z")

</div>

> [@bnicenboim](#):
>
> it would make more sense to have another package where tools are implemented for both interfaces

Yeah I thought about that right after I posted! I think that’s probably a better way to go. We really do need an SBC package. Would you be interested in collaborating with us on that?

---

<div class="post-metadata">

**Author:** ![bnicenboim](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bnicenboim/32/21197_2.png) [@bnicenboim](https://discourse.mc-stan.org/u/bnicenboim)\
**Post date:** [October 27, 2020, 8:04pm UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/5 "2020-10-27T20:04:45Z")

</div>

As you might have noticed by the amount of issues I’ve opened, and my recent posts :), I’m trying to understand SBC and trying to validate a quite complex model.  
If someone is leading this new package, I can definitely collaborate.

---

<div class="post-metadata">

**Author:** ![bnicenboim](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bnicenboim/32/21197_2.png) [@bnicenboim](https://discourse.mc-stan.org/u/bnicenboim)\
**Post date:** [October 28, 2020, 2:20pm UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/6 "2020-10-28T14:20:27Z")

</div>

I keep updating the functions, I might as well upload it with a package template. Should I do it in my own github for the time being? Or do you want to create a SBC project in stan’s git?

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [October 29, 2020, 8:25pm UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/7 "2020-10-29T20:25:37Z")

</div>

Ok great!

> [@bnicenboim](#):
>
> Should I do it in my own github for the time being? Or do you want to create a SBC project in stan’s git?

To start with I think do it in your own GitHub repo and then we can move it to stan-dev once it’s mature (this is what we did with other packages like the posterior package, which started out in my own personal github and then moved to stan-dev).

---

<div class="post-metadata">

**Author:** ![torkar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/torkar/32/1355_2.png) [@torkar](https://discourse.mc-stan.org/u/torkar)\
**Post date:** [November 7, 2020, 3:08pm UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/8 "2020-11-07T15:08:01Z")

</div>

@bnicenboim, I for one already look forward to having an SBC package! :)

---

<div class="post-metadata">

**Author:** ![hyunji.moon](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/hyunji.moon/32/9781_2.png) [@hyunji.moon](https://discourse.mc-stan.org/u/hyunji.moon)\
**Post date:** [November 8, 2020, 6:11am UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/9 "2020-11-08T06:11:39Z")

</div>

I wish to collaborate on sbc package. I have cmdstanr SBC code in [here](https://github.com/hyunjimoon/SBC_sensitivity/blob/master/laplace-approx/script/tools/sbcTools.R).

Usecases in for the codes are [here](https://github.com/hyunjimoon/SBC_sensitivity/blob/master/laplace-approx/script/prostate_cancer_sbc.Rmd) with the following summary.

```
sbc_res_bern_lik <- sbc_new(mod, modelName, data = data, N = N, L = L, save_progress = delivDir)
ppc.sbc(sbc_res_bern_lik,modelName, data)
uniformity.sbc(sbc_res_bern_lik, modelName)
plot.sbc(sbc_res_bern_lik, modelName, data)
plot.sbc_ecdf(sbc_res_bern_lik, modelName)

```

Resulting plots could be found in this [paper](https://github.com/hyunjimoon/SBC_sensitivity/blob/master/laplace-approx/writeup/Diagnosing_embedded_Laplace_approximation.pdf). Two things to note.

1. Code for plotting ECDF is included `plot.sbc_ecdf`, but code for ECDF difference ( following plot) is needed.  
 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/e/e895b891f684e366ab1dd9a44b4314408c9f25a3.png)
2. I felt the need for diverse measures for SBC other than chisq so I made `uniformity.sbc` to return Wasserstein and KM distance on top of chi-square p values.

```
MW1 <- function(bin_count){
  bins <- length(bin_count)
  unif <- rep(1/bins, bins)
  M <- sum(bin_count)
  tempf <- Vectorize(function(i) abs(bin_count[i]/M - unif[i]))
  val <- integrate(tempf,1,bins, rel.tol=.Machine$double.eps^.05)$value
  return(val)
}
MKM <- function(bin_count){
  bins <- length(bin_count)
  diff <- abs(mean(bin_count) - bin_count)
  val <- diff[which.max(diff)] / mean(bin_count)
  return(val)
}
MChisq <- function(bin_count){
  return(chisq.test(bin_count)$p.value)
}

```

---

<div class="post-metadata">

**Author:** ![bnicenboim](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bnicenboim/32/21197_2.png) [@bnicenboim](https://discourse.mc-stan.org/u/bnicenboim)\
**Post date:** [November 8, 2020, 7:41am UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/10 "2020-11-08T07:41:35Z")

</div>

That’s really cool!, thanks!  
I’ll set the package up this week and I’ll announce it here.

---

<div class="post-metadata">

**Author:** ![hyunji.moon](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/hyunji.moon/32/9781_2.png) [@hyunji.moon](https://discourse.mc-stan.org/u/hyunji.moon)\
**Post date:** [November 8, 2020, 9:59am UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/11 "2020-11-08T09:59:20Z")

</div>

@bnicenboim What additional features do you think are needed other than `ECDF_diff` mentioned above? To be honest, I had no difficulty doing SBC with the current code. Do we need to set up a new package? I guess we would need @jonah's opinion as well.

Happy to meet SBC colleague!

---

<div class="post-metadata">

**Author:** ![bnicenboim](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bnicenboim/32/21197_2.png) [@bnicenboim](https://discourse.mc-stan.org/u/bnicenboim)\
**Post date:** [November 8, 2020, 10:30am UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/12 "2020-11-08T10:30:04Z")

</div>

The idea was to create a new package:

> [@jonah](#):
>
> Yeah I thought about that right after I posted! I think that’s probably a better way to go. We really do need an SBC package. Would you be interested in collaborating with us on that?

and to make it compatible with both rstan and cmdstanr. Re new features, I agree that it’s probably enough. I’m doing SBC with even less.

if you want to take over and set it up, sure go ahead!

---

<div class="post-metadata">

**Author:** ![hyunji.moon](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/hyunji.moon/32/9781_2.png) [@hyunji.moon](https://discourse.mc-stan.org/u/hyunji.moon)\
**Post date:** [November 8, 2020, 11:08am UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/13 "2020-11-08T11:08:28Z")

</div>

Oh, I was not aware of that. I thought just Rcode would be enough like the current rstan.  
To provide samples for both interfaces, which language should be used? Is it correct that we cannot use rstan nor cmdstanr? Or could we store results obtained from one language and transfer to other language via functions like `$save_object`?

> [@bnicenboim](#):
>
> If someone is leading this new package, I can definitely collaborate.

@bnicenboim If I set up the package would you be willing to collaborate? For example, I have not followed with your [updates in rstan](https://github.com/stan-dev/rstan/commit/4bea50af1816af73d70fe07bcceaff9bc989a246) and parallelization.

---

<div class="post-metadata">

**Author:** ![bnicenboim](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bnicenboim/32/21197_2.png) [@bnicenboim](https://discourse.mc-stan.org/u/bnicenboim)\
**Post date:** [November 8, 2020, 12:51pm UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/14 "2020-11-08T12:51:42Z")

</div>

I don’t really follow. The functions should be able to take objects from rstan or cmdstanr and make the same objects of the class sbc.

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [November 9, 2020, 10:08pm UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/15 "2020-11-09T22:08:36Z")

</div>

@hyunji.moon Glad you’re interested in this too!

I can see a few options here. The first is an sbc package in the style of the loo package, where other packages provide their own methods for interfacing with loo. For example rstan has a loo method for stanfit objects, rstanarm and brms have loo methods for stanreg and brmsfit objects, etc. The loo package doesn’t know anything about these other packages, and relies on the other packages to provide it what it needs to do the computation.

Another option is for an sbc package that natively supports stanfit objects from rstan and CmdStanMCMC objects from cmdstanr. If the sbc package knows how to manipulate those objects internally then it could do what it needs to do without requiring the other packages to implement methods.

Regardless of the approach, we can make use of the existing code that’s in rstan (or rewrite it) and also the code that you all have written for cmdstanr.

What do you all think?

---

<div class="post-metadata">

**Author:** ![bnicenboim](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bnicenboim/32/21197_2.png) [@bnicenboim](https://discourse.mc-stan.org/u/bnicenboim)\
**Post date:** [November 10, 2020, 7:27am UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/16 "2020-11-10T07:27:52Z")

</div>

After looking at the sbc code for stan, I think that ideally an `sbc` function should just take a matrix/array/dataframe of generated data (the parameters\_ in the actual implementation) and a matrix/array/dataframe of posteriors and do the rank in R rather than in Stan. Then the function can work for rstan, cmdstanr, JAGS, etc, and one can even generate the data in R. (Now for example it’s easier to fit a truncated distribution in Stan than to generate it). Also we can avoid using regex that IMHO can miss cases (as it happened with my [commit](https://github.com/stan-dev/rstan/commit/fd704e473b220ba6e5d834c5b35ceb293f7980b7)).

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [November 10, 2020, 6:43pm UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/17 "2020-11-10T18:43:27Z")

</div>

Sounds good to me!

---

<div class="post-metadata">

**Author:** ![bnicenboim](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bnicenboim/32/21197_2.png) [@bnicenboim](https://discourse.mc-stan.org/u/bnicenboim)\
**Post date:** [November 10, 2020, 7:15pm UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/18 "2020-11-10T19:15:26Z")

</div>

@hyunji.moon, should I give it a try?

---

<div class="post-metadata">

**Author:** ![hyunji.moon](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/hyunji.moon/32/9781_2.png) [@hyunji.moon](https://discourse.mc-stan.org/u/hyunji.moon)\
**Post date:** [November 10, 2020, 9:45pm UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/19 "2020-11-10T21:45:12Z")

</div>

Sorry for the delay. Would it be ok if I let you know by today? I wish to know what the requirements are and looking into other packages, loo, [projpred](https://github.com/stan-dev/projpred), [posteriordb](https://github.com/MansMeg/posteriordb) for example.

---

<div class="post-metadata">

**Author:** ![hyunji.moon](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/hyunji.moon/32/9781_2.png) [@hyunji.moon](https://discourse.mc-stan.org/u/hyunji.moon)\
**Post date:** [November 11, 2020, 2:29pm UTC](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931/20 "2020-11-11T14:29:53Z")

</div>

@Dashadower and I have started SBC in this [repo](https://github.com/hyunjimoon/SBC). We wish to discuss the following:

General functions that the package could provide is just what @bnicenboim said, calculating the ranks and drawing the plot.

Sampling from the prior and refitting simulated values still need to be manually done by the user ([this](https://github.com/hyunjimoon/SBC/blob/bd5addd3464711f0a2d8d466ea3815eb47562c46/sbc_array.R#L34) process from our testcode), unless we parse the stan code (which could be unreliable). The first approach from @Bob_Carpenter’s [answer](https://stackoverflow.com/a/57718755) could be useful.

[Next page](https://discourse.mc-stan.org/t/sbc-for-cmdstanr/18931.md?page=2)