Multimodality reported in hierarchical prey-predator

[hyunjimoon]

Previous attempts to fit hierarchical prey-predator model in stan discourse suggests multimodality, which can be detected by n_eff <= n_chains. @tomfid

[hyunjimoon]

Could this be replicated in Vensim?

[hyunjimoon]

Michael suggests incremental precision as a remedy here;

You can try increasing the absolute and relative precision of the integration, but you’re fighting exponential errors so it’s a losing game. Chaotic systems are hard.

[tomfid]

I'm having a hard time making sense of what they've done here - multiple independent populations, or ... ? There's a lot going on in that thread.

Chaos and ODE precision issues seem like a distraction over realistic time horizons, unless maybe you're doing microbiology or something. Also, if the populations are independent, it's not chaotic.

[tomfid]

In the 4+2+1 parameter version in Vensim, there's no sign of multimodality - yet (but also no hierarchy yet).

This is MCMC:

And this is endpoints of random-restart Powell optimizations:

There's a surprising amount of dispersion in the Powell result - I guess finer tolerances would be in order. But since the startpoints are fully dispersed around the parameter hypercube, and everything ends up in roughly the same place, there's no sign of trouble here.

The alpha/beta/gamma/delta parameters do have substantial correlations.

[tomfid]

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.

[tomfid]

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.

[hyunjimoon]

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?

[tomfid]

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?

[tomfid]

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.

[hyunjimoon]

With the following settings:

N = 20
precision ={
    "S": 3, # # of draws from prior
    "M": 10, # # of draws from posterior (# of chains * # of draws from each chain)
    "N": N, # # of observation
    "Q": 3, # # of target_simulated_stock
    "R": 2, # # of subgroups for hierarchical Bayes
    "time_step": .03,
    "integration_times": np.arange(0, N) * .03 + 0.01,  # initial_time < integration_times, .03 need not be time_step, but length should be N
}

Ten minutes have passed, but there is no progress in HMC. Even with two estimated parameters, posterior geometry must be pretty bad.

07:33:39 - cmdstanpy - INFO - compiling stan file /Users/hyunjimoon/Dropbox/stanify/stan_files/pp_hierORsbc/pp_hierORsbc_data2draws.stan to exe file /Users/hyunjimoon/Dropbox/stanify/stan_files/pp_hierORsbc/pp_hierORsbc_data2draws
07:33:54 - cmdstanpy - INFO - compiled model executable: /Users/hyunjimoon/Dropbox/stanify/stan_files/pp_hierORsbc/pp_hierORsbc_data2draws
07:33:54 - cmdstanpy - INFO - CmdStan start processing
chain 1 |          | 00:00 Status
chain 2 |          | 00:00 Status

For comparison, it takes 30 seconds for two chains with R = 1 (no hier.)

chain 1 |          | 00:00 Status
chain 2 |          | 00:00 Status
chain 1 |█▋        | 00:01 Iteration:    1 / 1005 [  0%]  (Warmup)
chain 1 |██▌       | 00:06 Iteration:  100 / 1005 [  9%]  (Warmup)
chain 1 |███▎      | 00:09 Iteration:  200 / 1005 [ 19%]  (Warmup)
chain 1 |████▏     | 00:14 Iteration:  300 / 1005 [ 29%]  (Warmup)
chain 1 |█████     | 00:17 Iteration:  400 / 1005 [ 39%]  (Warmup)
chain 2 |█████▊    | 00:18 Iteration:  500 / 1005 [ 49%]  (Warmup)
chain 1 |█████▊    | 00:20 Iteration:  500 / 1005 [ 49%]  (Warmup)
chain 1 |██████▋   | 00:23 Iteration:  600 / 1005 [ 59%]  (Warmup)
chain 1 |███████▌  | 00:25 Iteration:  700 / 1005 [ 69%]  (Warmup)
chain 2 |█████████▏| 00:26 Iteration:  900 / 1005 [ 89%]  (Warmup)
chain 1 |██████████| 00:29 Sampling completed
chain 2 |██████████| 00:29 Sampling completed    

[hyunjimoon]

Tracing prior tightness, hierarchy, multimodality

Every run takes around 5 hours. Multimodality (if this is the cause) with lower level of delay will tame the posterior geometry, so I will try with higher value for birth fraction (through prior).

The first is one parameter estimation, with larger and tighter value for prior on fraction (lower dealy), with no hierarchy. Much faster than the other two with estimated parameter and hierarchy added (five subscript).

The problem is, larger prior mean destroyed the oscillatory reference mode I've been experimenting. I really should move off from this volatile model... and start inventory with one parameter estimation.

[tomfid]

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

[hyunjimoon]

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

[tomfid]

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.

[hyunjimoon]

The cause of long estiamtion time was, process noise wasn't off during estimation. When I added draws2data_s['process_noise_scale'] = 0.0, it is orders of magnitude faster.

[hyunjimoon]

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

[hyunjimoon]

After updating to process noise off during estimation, it takes 3 hours (no hierarchy and dot01 mpnoise).

I also fitted the hierarchical model and its information is stored in prey_predator_S1N200Q2R3/sbc_mpdot01.nc. Out of four chains, one chain ended up divergency, and the posterior space it explored needs tracing. An inventory model with four chains takes 10 hours and 70% is completed on average.

[hyunjimoon]

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.

[hyunjimoon]

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?

[tomfid]

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.

[tomfid]

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.

[hyunjimoon]

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?