Vensim -> InferenceData/Arviz

[tomfid]

Just FYI, I've started working on this, after improving a few minor things in the Vensim MCMC implementation.

Looking at the docs, it seems like this should be primarily on the Python side. It think Vensim should generate a couple output files that make import into InferenceData easy, probably via Pandas dataframes (because that's an easy way to import tab or csv).

As a sanity check:

The InferenceData schema consists of a number of 'groups' each of which is an n-dimensional table, which may have a couple of standard dimensions ('chain', 'draw', 'sample') as well as some model-specific dimensions. Not all of these groups will be relevant for Vensim.

Posterior is easy - this is basically the runname_MCMC_sample file we now create.

unconstrained_posterior - not relevant

sample_stats - seems to be mainly HMC stuff - not relevant - a few things like acceptance rate correspond, but not sure it's worth the trouble initially.

log_likelihood corresponding with the posterior samples - also easy. Currently this is in the _points file, but that format isn't really convenient, so it would be better to create a separate file or column in the _sample file to make import easier.

log_prior corresponding with the posterior samples - also easy. However, for this and the previous item, we'll need to introduce some separation in the payoff (.vpd) so that the prior and data likelihood can be distinguished. Currently we don't identify which is which. This is mainly a UI question.

posterior_predictive - this will take a little automation, to generate a sensitivity run using the posterior sample. Should be pretty easy. I think existing savelist functionality can be used to identify the relevant variables to include (for large models, definitely not everything). Seems like the mindset here is posterior predictives corresponding with the observed_data, but policy variables would also be interesting.

observed_data - I think this might be messy, because the data could enter the payoff directly via a *C item, or the user might have some custom structure in the model for model-data comparisons. In the latter case, it wouldn't be easy to automatically identify the data and include it. Also, Vensim data is typically going to have nonuniform time axes and big variations in dimensionality. This seems like a version 2.0 item.

constant_data - Also messy, for the same reasons. Would be good to have for documentation and comparison purposes, but storing everything is redundant, and could run towards a gigabyte or more for a big model. Maybe have a UI for storing only policy relevant changes that reflect experimentation, not defaults?

prior - a prerequisite for doing this is identification of prior components from the payoff (as in log_prior above). Not sure Vensim has a strong motivation for generating samples from the prior by themselves.

prior_predictive - combines issues of the posterior_predictive and prior above.

predictions - out of sample predictions seem like they have definite utility. For Vensim, normally this is going to mean "predictions about the future" and it seems a little weird to have them in a separate table when really they are just extensions of the time dimension.

predictions_constant_data - Like constant_data, good for documentation but maybe not essential right away, since the model/vdfx also document this.

[hyunjimoon]

thanks for sharing. i'm so glad. agenda2 for today would be to brainstorm a heuristic rule for issue vs discussions in this repository. we can use action trigger so that opening / closing the issue can automatically update the project schedule. charlie asked to have suggest some milestone for the tools i'm building and it'd be really helpful to chat about potential goal / milestone of this project.

relevant threads are

• 🎶 Benefits of Connecting Vensim to Stan #21
• Recommend using Xarray/InferenceData #46

i added this repository to GitHub Copilot Chat and am experimenting what it does.

[tomfid]

Appreciate the links - helpful stuff.

My heuristic would be to create an issue when there's some actual code proposed, particularly the python side for injecting Vensim output into InferenceData. For Vensim changes, I'll probably create a design ticket in our internal system (Assembla). So far I feel like I'm still just circling the problem.

[tomfid]

I'm still pursuing the approach of making this mostly Python side, but most of the work has to start on the Vensim side. So far:

• I've added a syntax for specifying priors to the optimization control file (.voc), so now you can write something like:
  0<=fraction alloc to process<=1, BETA(5,5)
  or
  0<= transmission rate[region] <= 5, LOGNORMAL(1,regional variation SD)
  0<= regional variation SD <= 10, CAUCHY(0,1)

• There's a new type of optimization payoff element for priors, so you can specify a prior in the .voc (for a parameter with a common distribution) or in the .vpd (for something complex, where you want to write an equation)

• The payoff now computes the prior terms and data likelihood separately. This enables sampling from the prior (not implemented yet, but easy.)

• Calibration payoff elements (i.e. terms in the data likelihood) can now generate synthetic data, drawing from the assumed error distribution. This is really the key piece enabling something like SBC.

• The MCMC algorithm has been upgraded to DREAM-ZS (still no gradients, but seems to be a big improvement - at least 2x convergence speed)

• Sensitivity runs can now ignore some reserved words, so we can make the sample more informative (for InferenceData) without screwing up reuse of the sample for projections.

[tomfid]

There's a little more to do on the Vensim side:

• Interface for priors
• Sampling from the prior, both for experimentation and synthetic data generation, and for MCMC initialization
• Output file format changes supporting InferenceData (so it's easy to import the Vensim output into a dataframe and convert that into InferenceData)

[tomfid]

There's at least one imminent question: the InferenceData schema includes a number of 'groups' that are tables for particular kinds of output. If you're importing from Stan or some other known source, things are directed to the right group. However, when importing from dataframes from an unknown source, the documentation/examples are extremely vague about how you direct the incoming data to the right group.

[hyunjimoon]

perhaps of interest on automated sbc idea

4. SBC_Who

User-Program WF: One actor and its tool

• knowledge to distribution
• Data and parameter parsimony paysoff
  User's main role is to translate managerial prior knowlege to statistical prior distirbution. Program's role is to translate translate explicit in

Step	Goal	Machine, work, command (P-rows have .function(input))	Human's work	Main Bayes Component	output format
1	Specify Deman_prior		H1.a Listing policy functions		.vpd
2	Specify Structural_prior	Vensim assists H1.a()	H2a. Translate mental model to SD model (behavioral classification)	Prior (parameter)	.mdl
3	Specify Variational_prior for parameter	PySD assists H2.a()	H3a. Classify parameters to est_param and ass_param	Prior (parameter)	.json
			H3b. Supply value for assmed_param
			H3c. Supply value for assumed_param_ts
4	Specify Variational_prior for data		H4a. Classify states to obs_state Select obs_state among stocks	Prior (data)
5	translate Manag. to Stats. lang.	M1.a PySD, .build_function_block(H2.a)			relation.stan
6	Specify_explicit		H5a. Choose family(:= dist. of y, penalty distribution for error)	Likelihood	draws2data.stan gq block, data2draws.stan model, gq block
			H5b. Choose prior_family (default: PERT Normal)
7	Specify_implicit		H5a. Supply {min, mode, max} value for est_param's prior param		draws2data.stan gq block, data2draws.stan model, gq block
			H5b. Supply sign (real, non-neg) for est_param's prior param
			H5c. Supply type (disc/cont) for est_param's prior param
8	translate prior knowledge to distribution		M2a. map U5.abc to distributions (PERT normal, Poisson)		$\theta \sim Normal(3, 1.5^2), \sigma \sim Gamma(5,2)$
9	predict	M3.a draws2data.stan, fit_prior_data.sample(), fit_prior_data = (U2.ab, U3.ab, U4.ab): Prior predictive check (opt-out prior)
10	infer to verify	P4.a Stan, data2draws.stan,.create_stan_program(U1.a, U2.abc, U3.a, U4.ab, U5.abc): Infer parameter from (synthetic) data (Test or autoCalib)			Prior predictive check plot (summary stats.)
11	Specify_tolerance		H6a. Set precision with iter_sampling (:= # of samples)		$\gamma$ from SBC-graphics
			H6b. Select posterior approximator from [[diagram(universal_adversarial)]]	Posterior_approximator
12	infer to validate	P5.a Stan, fit_post_draws.sample(), fit_post_draws = (P1, U3.ab, U4.ab, U5.ab): Posterior predictive check (opt-in prior)			Posterior predictive check plot

Bayesian WF: Between three modules

• Bayesian workflow paper
• more focused on SBC as prior calibrator; relevent references on theory and application can be found here

SilkRoad WF: Between six softwares (Vensim-SDA, Stan-SBC, BATS-SOPS)

• Program perceived Demand with Vensim-SDA, Compute scientific Draws with Stan-SBC, Supply rationing Data with BATS-SOPS
• Diagrams and details for each Software is detailed here
• SilkRoad project aims to connect the above whose motivation is detailed here.

Step	Output	Software	Symbol	Description	opensource? (language)
1. Program Perceived Demand	a. Perceived Demand	Vensim	👁 Eye	Reads mental model, Translates to cyclic directed graph generator	X (has free version)
1.	b. Analyzed Demand	SDA	🧠 Brain	Finds dominant cycle of generator, Maps with system behavior	O (Mathematica, R)
-
2. Compute Scientific Draws	a. Computed Draws	Stan	🐴 Workhorse	Builds posterior space, Runs HMC, ADVI, BFGS for representative draws	O (Stan connected to Python, R, Julia)
2.	b. Verified and Validated Draws	SBC	👌 Test	Diagnoses graphically, Calibrates architecture, policy, parameter prior discriminator	O (R)
-
3. Collect Supporting Data	a. Theoretical Policy Parameter	BATS (communicating with Yaman, Gönenç)	🦇 Explore	Specifies policy parameter for demanded behavior	O (Python)
3.	b. Empirical Policy Parameter	SOPS (communicating with Erling)	🚀 One-shot	Optimizes policy one-shot in stochastic dynamic system	O (Powersim)
-
Iterate and communicate	Hub pysd, readsdr	PySD, readsdr	🗣 Language		O (Python, R)

[hyunjimoon]

fuller detail in SBC_dictionary.md