I pushed a spatial hierarchical version. https://github.com/Data4DM/stanify/tree/main/vensim_models/mng_pop/gen-est%206param%20spatial%20hier

This adds 4 regions with variety in the alpha/beta/gamma/delta parameters, loosely structured after https://betanalpha.github.io/assets/case_studies/hierarchical_modeling.html that @hyunjimoon linked. Part of the intent was to create the opportunity for funnels.

I also added diffusion among the regions. This is potentially pretty interesting, because it means there are two ways to get similar behavior among the 4 regions:

1. low variety in the vital rate parameters, and similar initial conditions, so the regions are on similar but independent trajectories
2. significant migration among the regions - it doesn't take much to couple the regional trajectories fairly strongly
   Hopefully this creates an added opportunity for multimodality, where the variety and migration parameters are confounded.

I didn't get around to including the migration parameter in the statistics yet, so any multimodality at present is a consequence of something else.

I haven't run it much yet, but I think it does show evidence of multimodality, or of pathologies that cause Powell to get stuck occasionally. I haven't tested MCMC either, but will get there soon.

Both of the plots above use Weka - https://waikato.github.io/weka-site/index.html There's nothing special about it - it's just a quick way to generate a scatterplot matrix (though a little annoying because it wants csv not tab input). There are better options in both R and Python, but I was in a hurry and didn't bother to work out the details.

The first one is a plot of the runname_MCMC_sample.tab file generated from an MCMC run in Vensim. It would actually be nicer to use the _points file, because it includes the payoff, but you'd have to filter by status code to exclude or color the rejected points.

The second one is the runname_endpoint.tab file generated from the Powell optimizer when random multistart is active.

I've always wondered the role of startpoint.tab, endpoint.tab. I thought they are produced only for recording purpose i.e. where to start again but you using them as similar to MCMC samples?

I also added diffusion among the regions. This is potentially pretty interesting, because it means there are two ways to get similar behavior among the 4 regions:

Now that I think about it, this doesn't necessarily create multimodality. It might just be a diagonal ellipse.

Also, in my version a lot of the parameters are lognormal not normal, which creates zero-avoidance. That might mitigate any funnel degeneracy?

Actually ... it looks like there's no particular relationship between the migration rates and variety, but there is a pretty strong tradeoff between the predator and prey migration rates.

Elsewhere there are strongish alpha-beta and gamma-delta correlations, but everything else looks like a round fuzzy blob.

I don't have much intuition for this.

The high frequency oscillation in the left chart is troubling - looks like integration error. What's happening there?

How could integration error lead to oscillation? (as I thought it could only be measurement or process noise)

There's an example at https://metasd.com/2017/03/dynamics-of-the-last-twinkie/ - look for "however if you parameterize f() poorly". You can also do something like SMOOTH( STEP(1,10), 0.66 ) with TIME STEP=1 I think. However, these are discrete oscillations with Euler integration, and shouldn't happen with a higher-order auto method. Whatever is happening here looks like a discrete oscillation, but may be something weirder.

Variation of fitting time is large, as some end in 3min. For one parameter estimation, all four chains took four hours.

Could the issue here simply be that the big if-else statement in the process noise lookup table is expensive? By the end of the sim, each call is iterating over the full n data points before returning.

Good point, will test this with smaller datapoints, N.

With 200 to 20 datapoints, 4 hrs is shortened to 4mins. Thanks! Hmm is there a better way to address if else when N gets bigger?

Not sure what control structures you're allowed to use in a function for the solver, but a bisection search is probably the fastest general method - something like: https://medium.com/self-training-data-science-enthusiast/bisection-search-1d34bd2987ea

Normally it's implemented as a recursive function or while loop, but it could also be implemented by reordering the if-else statement. Something like:

if x < 2
   if x = 0
     y = 3.5
   else # x = 1
     y = 2.1
else
   if x = 2
     y = 2.7
   else # x = 3
     y = 4.4

You'd need roughly logn layers of nesting.

Actually, the simplest thing, if the x-axis is uniform, is something like

xmin = 0
xmax = 20
yvec = { vector of points }
n = 21

xint = (xmax-xmin)/(n-1)
index = int( (x-xmin)/xint ) # should also check min/max bounds

y = yvec[index] 

This produces stepwise output like Vensim look forward, but would be easy to modify to interpolate.

Implementing bisection search seems unavoidable, as if else statement for stan function applies not only to process noise, but also to the input data (e.g. customer order rate).

Have you experienced any speed issue regarding lookup function? Could order fulfillment rate lookup function be speedup using bisection as well?