Help with usage of MCMC
-
- Member
- Posts: 22
- Joined: Mon Jan 21, 2019 11:41 pm
- Vensim version: DSS
Help with usage of MCMC
Hello everyone, I hope you are doing fine and safe. I am working with a VENSIM model of the Monarch Butterfly and want to get the uncertainty of 6 estimated parms in the model using MCMC posterior sampling. This is the first time that I use VENSIM's MCMC and I am having trouble reaching convergence. I've tried modifying the algorithm settings in several ways and after more than 1 Million runs I am still not able to reach convergence. Any suggestions on how to improve my convergence would be highly appreciated. Thanks a lot!
Re: Help with usage of MCMC
I think the first question is probably whether your payoff is properly weighted. With Normal distributions, it should have an expected value roughly equal to the number of points.
Can you share your setup - if not the model and data, at least the .vpd and .voc files and a payoff report (.rep)?
Can you share your setup - if not the model and data, at least the .vpd and .voc files and a payoff report (.rep)?
/*
Advice to posters (it really helps us to help you)
http://www.ventanasystems.co.uk/forum/v ... f=2&t=4391
Blog: http://blog.metasd.com
Model library: http://models.metasd.com
Bookmarks: http://delicious.com/tomfid/SystemDynamics
*/
Advice to posters (it really helps us to help you)
http://www.ventanasystems.co.uk/forum/v ... f=2&t=4391
Blog: http://blog.metasd.com
Model library: http://models.metasd.com
Bookmarks: http://delicious.com/tomfid/SystemDynamics
*/
-
- Member
- Posts: 22
- Joined: Mon Jan 21, 2019 11:41 pm
- Vensim version: DSS
Re: Help with usage of MCMC
Hi Tom,
Thank you so much for your reply. I checked the reading that you suggested to me over the FB page (Optimization and the Banana of Death) which gave me the idea of trying to use simulated annealing. I just started the sims 30 mins ago but it seems that it is moving in the right direction. I think I can still use the samples from the posterior for sensitivity analysis if I used SA, right? (I don't see why not)
Here are my files, what do you think?
Thanks
Thank you so much for your reply. I checked the reading that you suggested to me over the FB page (Optimization and the Banana of Death) which gave me the idea of trying to use simulated annealing. I just started the sims 30 mins ago but it seems that it is moving in the right direction. I think I can still use the samples from the posterior for sensitivity analysis if I used SA, right? (I don't see why not)
Here are my files, what do you think?
Thanks
Re: Help with usage of MCMC
You can definitely anneal prior to generating the sample - you just want to be sure that the temperature is down to 1 before the end of the burnin period.
I don't see your attachments.
I don't see your attachments.
/*
Advice to posters (it really helps us to help you)
http://www.ventanasystems.co.uk/forum/v ... f=2&t=4391
Blog: http://blog.metasd.com
Model library: http://models.metasd.com
Bookmarks: http://delicious.com/tomfid/SystemDynamics
*/
Advice to posters (it really helps us to help you)
http://www.ventanasystems.co.uk/forum/v ... f=2&t=4391
Blog: http://blog.metasd.com
Model library: http://models.metasd.com
Bookmarks: http://delicious.com/tomfid/SystemDynamics
*/
-
- Member
- Posts: 22
- Joined: Mon Jan 21, 2019 11:41 pm
- Vensim version: DSS
Re: Help with usage of MCMC
Excelllent, I'll do that.
I am not sure why the files didn't get attached
Here they are again
Thanks!
I am not sure why the files didn't get attached
Here they are again
Thanks!
Re: Help with usage of MCMC
OK ... here's the issue. At present you have a lookup containing the data, used to generate real_report. The model samples the overwintering population on day 30 to generate "report". The payoff compares the two and computes a sum of squares. So far so good.
Problem #1: the data in real_report isn't a data variable, so it's defined at every point in the model. That means your payoff is comparing the model and data about 40000 times, instead of ~30 times for the actual annual data points. This makes the scale of your payoff about 1000x too big. The simple solution would be to get rid of the lookup and use a real data variable for the annual points. (Obviously you can keep the lookup if needed for other purposes.) Then the payoff will be evaluated only where data exists. Since you sample the overwintering population on day 30 of each year, you could position the data points at day 30, 395, 760, etc. so they fall at the corresponding point in the simulation.
Problem #2 is that the calibration payoff element is using "1" for the scale parameter. That's fine if you're just using Powell to find the best fit, but it doesn't work for MCMC, for which the likelihood has to be properly scaled. Fortunately this is pretty easy.
First, modify the payoff to use the Gaussian distribution instead of the Normal. This does 2 things:
- it makes the scale parameter the standard deviation of the data, which is easier to work with
- it includes a couple extra terms in the likelihood for proper scaling and penalization of the std dev, so you can estimate it if you want
Then you have several options:
- If you have an estimate of the std dev of the data, you can use that - probably this would be time varying, so you'd need another data variable for it.
- If you don't, you can guess. You might use 10% of the overwintering population, for example.
- You can make this a parameter in your model and estimate it, if you want to get fancy.
What I often do for a quick and dirty test is to use the Logarithmic transform option on the payoff element, and set the StdDev = 0.1. That says that the error is about 10% of the data, and usually gets things into a ballpark that makes sense right away.
If things are working right, the expected value of the Values term in the payoff (column E in the payoff report) will be roughly the number of data points. Currently it's about 10000 times that. This causes MCMC to fail, because the probability of accepting a point that should like within the confidence bounds is EXP(-10000) instead of the EXP(-1) that you would get with proper weighting.
The payoff is -374789
The following are the components of the payoff:
Type Component Contribution Percent Values Params Used Skipped
*C report|Real report/1 -374789 100 -374789 0 39418 -39418
Problem #1: the data in real_report isn't a data variable, so it's defined at every point in the model. That means your payoff is comparing the model and data about 40000 times, instead of ~30 times for the actual annual data points. This makes the scale of your payoff about 1000x too big. The simple solution would be to get rid of the lookup and use a real data variable for the annual points. (Obviously you can keep the lookup if needed for other purposes.) Then the payoff will be evaluated only where data exists. Since you sample the overwintering population on day 30 of each year, you could position the data points at day 30, 395, 760, etc. so they fall at the corresponding point in the simulation.
Problem #2 is that the calibration payoff element is using "1" for the scale parameter. That's fine if you're just using Powell to find the best fit, but it doesn't work for MCMC, for which the likelihood has to be properly scaled. Fortunately this is pretty easy.
First, modify the payoff to use the Gaussian distribution instead of the Normal. This does 2 things:
- it makes the scale parameter the standard deviation of the data, which is easier to work with
- it includes a couple extra terms in the likelihood for proper scaling and penalization of the std dev, so you can estimate it if you want
Then you have several options:
- If you have an estimate of the std dev of the data, you can use that - probably this would be time varying, so you'd need another data variable for it.
- If you don't, you can guess. You might use 10% of the overwintering population, for example.
- You can make this a parameter in your model and estimate it, if you want to get fancy.
What I often do for a quick and dirty test is to use the Logarithmic transform option on the payoff element, and set the StdDev = 0.1. That says that the error is about 10% of the data, and usually gets things into a ballpark that makes sense right away.
If things are working right, the expected value of the Values term in the payoff (column E in the payoff report) will be roughly the number of data points. Currently it's about 10000 times that. This causes MCMC to fail, because the probability of accepting a point that should like within the confidence bounds is EXP(-10000) instead of the EXP(-1) that you would get with proper weighting.
The payoff is -374789
The following are the components of the payoff:
Type Component Contribution Percent Values Params Used Skipped
*C report|Real report/1 -374789 100 -374789 0 39418 -39418
/*
Advice to posters (it really helps us to help you)
http://www.ventanasystems.co.uk/forum/v ... f=2&t=4391
Blog: http://blog.metasd.com
Model library: http://models.metasd.com
Bookmarks: http://delicious.com/tomfid/SystemDynamics
*/
Advice to posters (it really helps us to help you)
http://www.ventanasystems.co.uk/forum/v ... f=2&t=4391
Blog: http://blog.metasd.com
Model library: http://models.metasd.com
Bookmarks: http://delicious.com/tomfid/SystemDynamics
*/
-
- Member
- Posts: 22
- Joined: Mon Jan 21, 2019 11:41 pm
- Vensim version: DSS
Re: Help with usage of MCMC
Thank you so much Tom. That makes a lot of sense!.
Hopefully, this will make everything run soft and smooth.
I hope you don't mind if I get back in touch with you if I come with any problems (hopefully not)
Thanks and stay safe!
Hopefully, this will make everything run soft and smooth.
I hope you don't mind if I get back in touch with you if I come with any problems (hopefully not)
Thanks and stay safe!
Re: Help with usage of MCMC
Please do - I'm curious to see how this works out.
/*
Advice to posters (it really helps us to help you)
http://www.ventanasystems.co.uk/forum/v ... f=2&t=4391
Blog: http://blog.metasd.com
Model library: http://models.metasd.com
Bookmarks: http://delicious.com/tomfid/SystemDynamics
*/
Advice to posters (it really helps us to help you)
http://www.ventanasystems.co.uk/forum/v ... f=2&t=4391
Blog: http://blog.metasd.com
Model library: http://models.metasd.com
Bookmarks: http://delicious.com/tomfid/SystemDynamics
*/
-
- Member
- Posts: 22
- Joined: Mon Jan 21, 2019 11:41 pm
- Vensim version: DSS
Re: Help with usage of MCMC
Hi Tom, I followed all your recommendations and things are starting to look considerably better. However, there is a very odd problem that I cannot wrap my head around and it is impeding me to get my posterior samples.
The thing is that the Acceptance Rate, Max PSRF, and CMCP are right on their ideal values almost from the very beginning while on the burning period. However, as soon as the burning period ends, the Max PSRF shoots through the roof and the Acceptance Rate keeps decreasing until reaching extremely low values. Why is that? and how can I fix it?
I thought that the burn-in period did not affect directly the walkers, just the sampling, right? or am I missing something?
Again, thanks a lot for your help
Best,
The thing is that the Acceptance Rate, Max PSRF, and CMCP are right on their ideal values almost from the very beginning while on the burning period. However, as soon as the burning period ends, the Max PSRF shoots through the roof and the Acceptance Rate keeps decreasing until reaching extremely low values. Why is that? and how can I fix it?
I thought that the burn-in period did not affect directly the walkers, just the sampling, right? or am I missing something?
Again, thanks a lot for your help
Best,
Re: Help with usage of MCMC
I think this has to do with a reset of the stats that happens at the end of the burnin period. Tentatively, the new values may not be right for a while. This is potentially ignorable, but I'd have to see a run to be sure. I guess I'd be most concerned about a drop in the acceptance rate. One thing to check is to be sure that the temperature is down to 1 at the end of the burnin period. (With the payoff reweighted, you may not need to anneal at all.)
/*
Advice to posters (it really helps us to help you)
http://www.ventanasystems.co.uk/forum/v ... f=2&t=4391
Blog: http://blog.metasd.com
Model library: http://models.metasd.com
Bookmarks: http://delicious.com/tomfid/SystemDynamics
*/
Advice to posters (it really helps us to help you)
http://www.ventanasystems.co.uk/forum/v ... f=2&t=4391
Blog: http://blog.metasd.com
Model library: http://models.metasd.com
Bookmarks: http://delicious.com/tomfid/SystemDynamics
*/
-
- Member
- Posts: 22
- Joined: Mon Jan 21, 2019 11:41 pm
- Vensim version: DSS
Re: Help with usage of MCMC
Hey Tom, I already ran sims with and without annealing, and with and without the burn-in period, but I still see the same; as soon as the samples begin to get drawn, the stats move from an almost-ideal position to an awful one.
I also made some runs with SA and making sure that temp goes down to 1 before sampling; still when the samples start to get drawn is when things go bad (see image)
I also made some runs with SA and making sure that temp goes down to 1 before sampling; still when the samples start to get drawn is when things go bad (see image)
Re: Help with usage of MCMC
It looks like your cooling period is longer than your burnin period, so there's an abrupt drop in temp from 5.77 to 1. I'd try changing the lengths so that doesn't happen, or just turn off the cooling altogether (you may find that you don't need it anymore).
/*
Advice to posters (it really helps us to help you)
http://www.ventanasystems.co.uk/forum/v ... f=2&t=4391
Blog: http://blog.metasd.com
Model library: http://models.metasd.com
Bookmarks: http://delicious.com/tomfid/SystemDynamics
*/
Advice to posters (it really helps us to help you)
http://www.ventanasystems.co.uk/forum/v ... f=2&t=4391
Blog: http://blog.metasd.com
Model library: http://models.metasd.com
Bookmarks: http://delicious.com/tomfid/SystemDynamics
*/