Getting stanify feedback from stats audience

[hyunjimoon]

Preparing for Stan discourse post and/or Andrew's blog about stanify

@jandraor and @hyunjimoon will develop three casestudies with @tomfid, @hazhirr as well as Stan community's leaders. It aims for Bayes workflow on prey-predator, SEIR, WIP-inventory management model.

Role models:

• Aki's birthday problem on gaussian process workflow here which builds eight different models depending on additional structure (long, short, event term) and distribution for variance parameter.

• Martin's Simulation-Based Calibration Checking for Bayesian Computation:The Choice of Test Quantities Shapes Sensitivity

• Michael's https://betanalpha.github.io/assets/case_studies/gaussian_processes.html#43_Approximating_Gaussian_Processes

The precision matrix can also be directly related to alternative representations of Gaussian processes, such as spectral densities and stochastic differential equations. Although these are more mathematically sophisticated representations they can offer a wealth of understanding about the functional behavior encoded within a Gaussian process, complementing the interpretability of the covariance function representation.

• Hazhir's Heterogeneity and Network Structure in the Dynamics of Diffusion (Agent based vs compartment model)

Comparing disaggregated DE models to AB models may be useful in designing stochastic compartment models that capture heterogeneity and network effects using the fewest additional compartments, if the network structure is known and stable. 3 Testing this proposal is beyond the scope of this paper.

• Angie's Predicting Engine Failure With Hierarchical Gaussian Process shows re-parameterization based on Stan's divergence plot

• King's Avoidable errors in the modelling of outbreaks of emerging pathogens, with special reference to Ebola, model excluding demographic noise biases estimation; CTMC and SDE

Format:

• for term consistency we will use the column of stan term defined in stanify which I plan to ask feedback from at least two community leaders from each side
• final format is quarto

Questions

1. effect of measurement noise and process noise
   a. relative scale between the two noise and their priors
   b. meaning and formulation of process noise in hierarchical context

2. examples that HMC can pass (fit fast) and/or fail with diagnostics (fail fast) i.e. joker combination 🃏 Joker that completes Batman #71 degeneracy from colinearity, hierarchy, multimodality and its remedy (Michael's managing degeneracy)

3. bad parameterization example detected by SBC: mixed order in gaussian mixture, standard-deviation vs precision formulation

4. types of hierarchy #66 demarcating the role of Vensim and Stan

5. hierarchical compartment model

• Niko's monster-model revisited
• addressing computational challenges with adaptive precision and time step

Goal

1. Leverage the authority of two community leaders during preparation (Hazhir and Tom, Andrew and Aki or Michael) to communicate the importance of the need for time and structure convergence of estimator / generator to 1. This stan post on hierarchical ode.

In Stan community, $\frac{t_{est}}{t_{gen}} > 1$, but in dynamic simulation community, $\frac{t_{est}}{t_{gen}} << 1$ where $t_{gen}, t_{est}$ are time spent in designing generator (parameter scaler values to state time-series matrix) and estimator (reverse of generator). From trustable model perspective, $\frac{t_{est}}{t_{gen}} \sim 1$ and symmetric generator and estimator are ideal. However, I discovered some model components that were asymmetrically focused in two communities. Topic of adaptable time_step, identifiablity, diagnostics are much valued in estimator-focused community where as process noise, aggregation get greater attention in generator-focused community. The difficulty of SBC lies in matching these asymmetry which were not surfaced when only one side of the joint model (joint of generator and estimator which forms joint distribution of parameter and observation) was the focus.

2. extend the boundary of model (add computation to managerial-statistics for SD people, add managerial to statistics-computation for Stan people) e.g. different perception of time-step asked in stan discourse

3. reparameterization strategies based on visual diagnostices for both decision-maker and computer

• e.g. from hierarchical gaussian process casestudy:

Three case studies

1. prey-predator

a. the importance of tight prior: sensitive parameter retrieval
-- replicating degeneracy observed in stanify (two birth frac. estimation, m_noise as beta(2,2), pnoise scale .01 in generator) in vensim

b. importance of tight measurement distribution (negative-binomial (mean-variance linear VS lognormal (too fat-tail; power of 2; mean is exp( µ + 2 1 σ 2), the variance is exp(2 µ ) exp(σ 2 )(exp(σ 2 ) − 1))

c. effect of process noise (process noise have oscillation-killing, aggregation effect which leads to stable but too much aggregated estimation which can be thought as layer 1 (question of new subscript can be discussed; new ship with its engine type included in train set VS new ship with its engine type not included in train set from the paper on korea navy engine failure)

d. clustered subgroup with normal mixture prior for population distribution

(Jair's comment: SDA is based on this model so changing dominant loop can be interesting)

2. SEIR

a. degeneracy from collinearity (mean generation time (1/gamma + 1/sigma), 2. basis reproduction number (gamma/beta))
b. interacting partial pooling (e.g. age group), ode function is called n_obs times
c. sampling distribution: negative binomial (poisson (overconfidence and biased estimates))

3. Inventory management

a. degeneracy from separating compartments degeneracy from multimodality (Hazhir discovered a bad parameter recovery with hierarchical Bayesian with powell)

b. non-interacting partial pooling documented in #17, ode function is called n_obs * n_region times

c. effect of timestep (decrease the time step to make sure it is not integration error, exponential decay going to 0; stock for expected customer order)

4. Future connection

• The role of stochasticity and heterogeneity with process noise #86 modeling capability of SDE in Vensim or Stan? Or do we have a walkaround (Jair says CTMC)? Michael's inferring SDE. This casestudy on hierarchical gaussian process.

• Ed's Combining stock-and-flow, agent-based, and social network methods to model team performance

the differences among the productivities of the diversity scenarios is not due to the number of members from each subgroup, which is constant among the three scenarios. Rather it is from members’ positions in the communications hierarchy.

[hyunjimoon]

How much quantitative data is needed to reliably infer parameters of an SD model?
to define three steps as follows:

1. parameterize SD model
   input: pre-processed differential equations (.mdl) • reparameterize to management friendly form: Tom's parameter taxonomy • reparameterize to statistical computation friendly form: express parameters using notations defined in Bayesian taxonomy
   output: structure (processed differential equations)

2. collect and construct driving data array
   input: structure + pre-processed data (data lake)

   conditional on structure
   • collect required data
   • frame data as data array (labeled multidimensional array)
   • frame prior knowledge on parameters as driving prior distribution
   output: consistent (structure, numeric)

3. flexible tests using experimental priors
   input: (structure, numeric) + managerial knowledge • choose test quantities: (log likelihood + log prior penalty) is default and modeler can add other test quantities • choose setting: different experimental priors for estimated, policy parameters
   output: numeric and graphical diagnostics
   ——
   Iterate 1,2,3 and inspect diagnostics
   • numeric: d(target_obs, target_simulated_obs), inferential and predictive empirical continuous rank probability score (defined in Michael’s writing), Rhat, effective sample size, p-loo, expected log predictive density for a new data point, E-BFMI • graphical: scattered divergence on pairs plot, sbc rank plotfinal output: useable consistent (structure, numeric, setting)
   ——stanify draws2data2draws automate the third step:

def draws2data2draws(structure, numeric, setting,  S, M, N):
  model = vensim2stan(structure)  model.set_numeric(numeric)

  model.set_setting(setting, N)
  prior_sample = sample(experimental prior, S)
  target_simulated_obs = draws2data(model, prior_sample)

  for s in range(S):
    posterior_sample[s] = data2draws(model, target_simulated_obs[s], M)
   
  def draws2data(model, prior_sample)
    return generate(model, prior_sample).target_simulated_obs
  def data2draws(model, target_simulated_obs, M)
    return estimate(model, target_simulated_obs, M)

  def diagnose(prior_sample, posterior_sample, target_simulated, target_simulated_obs, test_quantity):
    return compare(test_quantity(prior_sample, target_simulated, target_simulated_obs, target_obs),
                              test_quantity(posterior_sample, target_simulated, target_simulated_obs, target_obs))
  
  return diagnose(prior_sample, posterior_sample, target_simulated, target_simulated_obs, target_obs, ('loglik'))

——
I wish to get feedback on the overall procedure defined above and two new design decisions:

1. I removed the concept of assumed parameter as it can be framed as estimated with point mass prior
2. I separated driving prior (step 2) from experimental prior (step 3). From my experience, Stan users spend most of their time tweaking priors until they get good enough diagnostics. This computational consideration should be separated from our real knowledge of a parameter.

[hyunjimoon]

[hyunjimoon]

• Jair used reporting rate \rho \sim beta(2,2), \beta lognormal
• Angie used beta(2,2) for measurement noise prior which caused divergence; perhaps it was too big considering normal(.1, .01) didn't cause divergence.

[jandraor]

questions to ask to statistics experts via blog

What is the advantage of adding process noise to the SD prey-predator in comparison to the time-continuous Markov chain version?

[hyunjimoon]

Correlation endgoal

[hyunjimoon]

Both Stella and Vensim has output for xml.
i) python -> inferencedata (netcdf) -> R -> SBC
ii) vensim -> netcdf -> R -> SBC

Mike has a code that translates netcdf to rvars here. Explanation for rvars by Matthew and SBC seems to support rvars from here.

[hyunjimoon]

Thank you, Andrew Gelman and Tom Fiddaman @tomfid for supporting this project. Below is what @Dashadower and I planned together.

Hierarchical ODE SBC

Time range

12 weeks from Dec.W1 ~ FebW4.

Goal

1. Speed up
   measure of success:
   a. parameters combination (which pair and how many) for an inventory management model whose estimation ends in .5 hour and pass arviz diagnostics (rank, sbc, rhat, khat)
   b. at least one stan discourse question or asking feedback from ODE experts (Juho, Sebastian, Charles, Yi)

• trace bottleneck and experiment with what causes 13hours for inventory model estimation
• potential causes are less adaptive warmup (remedy for the monster model by Niko), diffuse prior, hierarchy, not using reduce_sum etc. Based on the auto-hierarchy stanfile translation and experiment setting Angie prepared, Shinyoung will experiment with different model settings to find the cause of the bottleneck and apply existing methodologies.
• .5 hours is based on the fitting time from Niko's report
• Multimodality reported in hierarchical prey-predator #89 (comment)

2. Connect
   measure of success:
   a. arviz function supporting sbc plot of a hierarchical model
   b. sbc ecdf_diff and coverage plot of applying a. to inventory model

• the decision to connect to R vs doing all in python is deferred to the end of December; this contains blueprint for connection to R
• regardless of the decision, Angie will develop inferencedata-based sbc plot for a hierarchical model by the end of Dec. ~ mid Jan. with @OriolAbril
• Shinyoung will apply this plot (or the logic - if we decide to connect with SBC package) to the inventory model after 1's completion.

3. Action
   measure of success
   a. replicate blog's output on the saved cost by 29 ("Expected value of including uncertainty = 29.31")
   b. add sbc, plot_ppc, pair plot and other basic diagnostics to blog's storyline(from 2)

• Angie will prepare vensim model by the end of Jan. Shinyoung and Angie will use this model to replicate blog after 1 and 2's completion.
• Bayesian decision flow on supply chain #105

Output

findings will be presented in quarto or notebook format which will be posted in Andrew's blog for feedback (as planned above).

[hyunjimoon]

12W4

1. difference between two pr? might be better for shinyoung to take the lead on refactoring; how long?
2. loglik plots? defining degeneracy diagnostics for dynamic model

• multimodal loglik compare: Degeneracy diagnostics for dynamic model #134

[hyunjimoon]

with @Dashadower next week

• sbc rank plot for data-loglik (M by S)
• sbc rank plot with lp__ and divergence mapping (hypothesis: the two have high negative correlation, meaning low lp__ would have create high divergence from bad geometry

[hyunjimoon]

Shinyoung, Angie 1W1

did

1. plot with xarray dataset (on time coordinate)
2. arviz assign within group and sum between prior_predictive and posterior_predictive

food for thoughts

1. variation ratio to each layer as sbc test quantity

With user-defined heterogeneity decomposition of total heterogeneity (variation) to individual (unit) level and population level e.g. 3:1 =log_pop: log_unit, can we do sbc with test quantity .25log_pop-log_unit?

2. could we relate loglike_unit and loglik_pop with contraction, containment, compromise?

i) contraction: unit_layer (likelihood) is more informative, but consistent with population_layer (prior)
loglik_unit << loglik_pop

ii) containment: population_layer (prior) is more informative, but consistent with unit_layer (likelihood)
loglik_unit >> loglik_pop

iii) compromise: tension btw information encoded in population and unit
loglik_unit ~ loglik_pop

3. example for heterogeneity contribution to multilayer (one layer added to unit-population above)

: manufacturer (hyper-population) + engine type (population) + engine (unit)
For forklift data we had from paper below,

Kim22_FailureFuncLogisticsBayesAgg.pdf

[hyunjimoon]

@Dashadower could we try using adjoint method? https://discourse.mc-stan.org/t/slow-hierarchical-ode-model/30001/2

[hyunjimoon]

1W2 Shinyoung, Angie agenda

1. merge shinyoung's branch with pdoc branch in pr
2. is it possible to translate vensim subscript?

[hyunjimoon]

Possible ideas to develop:

• Plot design for hierarchySBC stanify#6

[hyunjimoon]

2W4

• data (exo_demand) + tests
• loglik-based plot, not one theta as in https://gist.github.com/Dashadower/face0ff696d83a13d0ea1624b70c7917
• • meaning https://github.com/Data4DM/stanify/blob/68cde36549923ff29764171ab1dfc955bfb1bae9/stanify/builders/utilities.py#L255?

[Dashadower]

Data4DM/stanify#45

[Dashadower]

Data4DM/stanify#45

Especially supporting subscripted data/lookup