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.
For now parallelization is not working (https://github.com/stan-dev/cmdstanr/issues/326). But maybe this is useful for someone else.
Cool, thanks! We have an open issue for adding SBC functionality to CmdStanR
but haven’t gotten to it yet. Any interest in contributing this to the package?
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), 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).
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?
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.
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?
Ok great!
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).
I wish to collaborate on sbc package. I have cmdstanr SBC code in here.
Usecases in for the codes are here 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. Two things to note.
plot.sbc_ecdf, but code for ECDF difference ( following plot) is needed.
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)
}
That’s really cool!, thanks!
I’ll set the package up this week and I’ll announce it here.
@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!
The idea was to create a new package:
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!
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 I set up the package would you be willing to collaborate? For example, I have not followed with your updates in rstan and parallelization.
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.
@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?
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).
Sounds good to me!
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, posteriordb for example.
@Dashadower and I have started SBC in this repo. 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 process from our testcode), unless we parse the stan code (which could be unreliable). The first approach from @Bob_Carpenter’s answer could be useful.