How to use Markov Chain Monte Carlo in Vensim?

[gwr]

I am greatly interested in combining Bayesian Statistics and System Dynamics. Therefore I greatly appreciate the inclusion of MCMC within the release of Vensim 6.0 and would greatly encourage further development in that regard.

From reading through the instructions in the reference manual I am not quite sure how employing the MCMC-algorith in Vensim is done in order to approach Bayesian Statistics. Maybe you could help me therefore with these questions:

1. How can Prior Probability Distributions be included?

In Bayesian Statistics one usually has to assume some kind of Prior distribution for stationary (time invariant) parameters. These can be subjective or objective but in the end it means that one should be able to provide a prior probability distribution for any prarameter that is to be calibrated - e.g. for which a posterior distribution should be calculated given data [ prob( parameter | data, model ) ~ prob(data | parameter, model) x prob (parameter) ]. How can this be done in Vensim?

2. Why will MCMC in Vensim use optimization - e.g. what is the output of MCMC?

I am not quite sure, that I understand, why the MCMC - algorithm employed in Vensim uses the optimizer? As far as I know the full-fledged MCMC is a sampling algorithm that converges to sampling from the posterior distribution of the parameters [e.g. prob( parameter | data, model) ] after some burn-in-period. So the output in its general form should be a sample of points that are taken from the posterior-distribution of a parameter.

Optimization to my understanding of Bayesian methodology only enters the picture when we have to answer a decision like "which is the most likely parameter value?" (point estimation). Here, depending on a specific loss function, either the expected value of the distribution (Bayes estimate), the median value of the distribution or the maximum aposteriori estimate (MAP) can be used to answer that question. Optimization is used (to my knowledge) to arrive at the MAP-value which in the case of a constant prior value will be the maximum likelihood estimate. So could you elaborate here a bit more?

3. Are there any restrictions?

As for Kalman filtering and FIMLOF I am aware of the fact that it should not work to well with bigger models. Does MCMC have similar restrictions? What about Runge-Kutta-integration methods?

Kind regards,

Guido

[LAUJJL]

Hi Guido

I think too that Bayesian statistics is interesting as it permits to improve prior knowledge by mixing it with new experiences to deliver a new or posterior knowledge. It looks a bit like a negative adjustment loop. But I did not see any reference in the Vensim documentation about Bayesian statistics. The only time where Vensim is using some Bayesian methods is with calibration using likelihood functions.

From what I remember, Anylogic worked some years ago with a company specialized in Bayesian statistics and included in its package some Bayesian statistics. But it does not seem to have lasted very long and there does not seem to be any more references in Bayesian statistics.

There is a software that uses prior and posterior functions and works either in a static or dynamic way. This is Analytica. It proves that these functions could be implemented in Vensim.

Using Bayesian methods could help build a bridge between statistical methods based on data and SD based on causality.

Besides that I cannot answer to your questions and am too interested with any answer.

Best regards.
JJ

[gwr]

Hi JJ,

here it is to be found in the release note to version 6.0:

http://www.vensim.com/documentation/mcmc_sa.htm

It says: "The equilibrium distribution of the random walk reflects the probability density of the payoff. This makes it possible to directly assess the posterior probability of the calibration given the data and priors (if used)..." which suggests that using prior distributions should be possible; and after all, Markov Chain Monte Carlo is at the heart of Bayesian Statistical methods, isn't it?

So far I have experimented with Rudolf Kulhavý 's approach -- he nicely refers to it as a grey-box approach -- and rebuilt it using Mathematica (his source code is in R). He uses the particle filter which to my knowledge is an approximate Bayesian approach using importance sampling ("Bayesian Bootstrapping"):

http://www.systemdynamics.org/conferenc ... LHA346.pdf

As far as I understand this appraoch is not limited to linear dynamical systems as the Kalman-Bucy-Filter and might make a very powerful addition to the SD-toolbox. So I am very curious about the answer we will be getting here.

Best,

Guido

[LAUJJL]

Hi Guido

Thank you for the information. These features being not in the interface, it is easy to miss them.

All these sophisticated mathematic and statistical tools need to be tested not from an academic point of view, but from a practitionner one.

I have already difficulties to use properly the common Vensim tools and I am very prudent using these tools and would like to see them practically used.

Regards.

JJ

[gwr]

Hi JJ,

what would be the practioner's criteria in judging advanced (and very general) statistical tools in your view?

Following Kurt Lewin to me (at least in applied sciences) academia (=theory) and application (=practice) are not so much opposites but two sides of the same coin: Without theory what will you gain by experimenting with a model or with reality? (Just compare the speed of technological progress in history before mankind came up with the scientific method and after) So what is needed imo are ways to make the scientific method more easily applicable in day to day (or at least quarter to quarter) decisions.

Applying Bayesian statistics one would need to know a bit what one is doing otherwise you again cannot interpret what you are seeing correctly. On the other hand having practical use in mind I do find the concept of Bayesian statistics very appealing. You can map pretty much any relevant uncertainty be it measurement error or left out structure in equations or parameter uncertainty. Following Jaynes and others Bayesian prior probabilities would also not need to be simply subjective there are objective priors out there - but they nevertheless could be subjective as well allowing for a practioner to model his "gut feelings".

Whenever I do show sensitivity runs that is what people find most appealing. So then I find that having sensitivity runs where your probability distributions are tested and consistently modified according to the data you have at hand a quite powerful concept. How else to do proper risk modeling in a consistent way not prone to cognitive distortions? Following what Kulhavý lays out I cannot see anything but very practical utility of these features.

Regards,

Guido

[gwr]

Being quite "impatient" I found this posted by Tom Fiddaman on the Vensim website:

http://vensim.com/wp-content/uploads/20 ... -C-web.pdf

Now, while this is not really a substitute for some more elaborate answer, I do get some of the picture from what is presented under the topic "extended payoff". Obviously one can hav a penalty function for the payoff which is of course how one could interpret a prior distribution function which is multiplied by the likelihood in Bayesian inference. This is further explained in the documentation:

http://www.vensim.com/documentation/ext ... ayoffs.htm

If I understand this correctly, then prior information on parameters might be given by an inequality constraint that is computed at the initial time only. That would mean that one might only define uniformly distributed priors (e.g. where my prior knowledge is that the parameter is likely to be in a finite range but I do not have any further information on the distribution within that range) - or am I wrong here? What if I wantede a beta distribution or a truncated normal distribution or anything else for that matter?

Again, more information desperately needed.

Guido

[LAUJJL]

Hi Guido

I am very suspicious with using too sophisticated mathematical tools. I am mainly a business man, although having studied extensively theoretical mathematics and applied statistics in the University prior to studying management. I always try to be practical and beware of academic theories. I like too to understand my models like 1 + 1 = 2 and to much sophistication especially black box one's, will restrict that understanding.

I do not know the kind of studies you are on, but it seems to me that before wanting an answer one must first consider the problem studied. One idea may be good in one situation and bad in another.
I personally prefer to study first the problem and after to eventually study some new or old or forgotten tools.

To resume, for me there is no universal truth at least in social sciences and especially in business.
Any preconceived idea or method will limit my freedom of thinking.

So I cannot answer your questions unless you give me a concrete example.

Best regards.
JJ

[gwr]

Hi JJ,

sorry to have caused misunderstanding, but I have put forward a question to be mainly answered by the Vensim support, as it concerns technical issues that must be answered by people that have implemented them. I might be able to give you an example for Bayesian reasoning but that needs answering the technical questions first. What do you find misleading or impractical in the approach presented by Kulhavý ? (Not the mathematical formulas but his description of the benefits)

I do value business but I do also value "truth". That may be some old fashioned idea from philosophy but in the end I know in how many ways people can err due to cognitive biases and distortions. I still do not know any other method that might reliably provide truth than the scientific method. What should I trust, revelations or "intuition" instead? Following Mr. Kahneman just a bit I would be careful. Isn't that one reason we start modeling in the first place: To make our assumptions explicit and somehow testable?

If I remember correctly, than it was you, who asked about optimizing large models with ´000s of parameters. Now, I personally would be much more careful about optimization than I would be and am about applying Bayesian Statistics in consistently accounting for uncertainty. A matter of taste maybe?

Let us wait what Tom or Tony have to say to the questions posed.

Kind regards,

Guido

[tomfid]

I'm traveling, but hope to say a bit more about this soon.

[gwr]

I'm traveling, but hope to say a bit more about this soon.

I am not giving up hope either.

Guido

[tomfid]

Sorry for the delay. I've been holding off to compose a better answer and to chase a small bias in some results, but perhaps I shouldn't let the perfect be the enemy of the good.

I'll clean up and post some tutorial materials from the SD conference and my test cases as soon as I'm back from Boston this weekend. There are actually 2 topics here: MCMC as an algorithm, and Bayesian approach to estimation.

On a related note, there's a new random number generator coming in the next release.

[tomfid]

Coming back to the original questions:

1. Priors

The mechanism for including priors is to include them in the payoff.

This has actually always been possible in Vensim by creating a policy payoff that computes a likelihood, possibly combining likelihood of data points and priors. Something like:

Code: Select all

payoff = prior likelihood + point likelihood
prior likelihood = IF THEN ELSE( Time = INITIAL TIME, <expression for parameter likelihood>, 0)
point likelihood = IF THEN ELSE( data var <> :NA:, <expression for data point likelihood - e.g. -((model var - data var)/err sd)^2>, 0)

with payoff file

Code: Select all

*P
payoff/1

This gets a bit easier in many cases with the new extended payoffs, because you don't have to write equations - you can implement much of the above in the payoff definition.
http://www.vensim.com/documentation/ext ... ayoffs.htm
You might have something like:

Code: Select all

prior likelihood = IF THEN ELSE( Time = INITIAL TIME, <expression for parameter likelihood>, 0)

with payoff file

Code: Select all

*C
model var|data var/weight
*PF
prior likelihood/1

2. MCMC

The output of MCMC is the sample of points visited by the random walk. Assuming good behavior, the properties of the sample should represent the posterior distribution. In many cases, the data will be the dominant contributor to the posterior, so whether you define priors or not is irrelevant, but that's not always the case.

MCMC isn't strictly an optimizer, but it uses the same infrastructure as the hill-climbing optimizers - you specify a payoff and a list of parameters, plus perhaps some control settings. Also, by adding an annealing schedule, the algorithm does simulated annealing, in which case you might ignore the distributional results and only pay attention to the maximum.
http://www.vensim.com/documentation/mcmc_sa.htm

3. Restrictions

There are no restrictions on the model, and no additional computational burden (unlike Kalman filtering). However, it doesn't do anything about state estimation & updates, so there's still a role for the Kalman filter.

[gwr]

Tom,

thank you very much for your answer which is indeed helpful. since I am mostly interested in employing the Bayesian approach to model calibration I will have to set up a payoff that will combine a prior and a likelihood as a product if I remember this correctly so the sum given here might only work if I take logarithms on both sides?

I will hopefully find some time to work out the example given by Rudolf Kulhavý using Vensim's MCMC as there should be comparable results. As to state updates this is also addressed in the Kulhavý paper as driving and measuring noise are considered. Maybe my understanding here is rudimentary but Particle Filtering (the appraoch taken by Kulhavy) is to my knowledge an approximate Bayesian procedure (e.g. importance sampling) but maybe this topic belongs to the System Dynamics forum where there are plenty of academics to fill in the blanks.

Kind regards,

Guido

[tomfid]

Yes, you're right, I should have said log likelihood throughout.

I'll be interested to hear how things work out. I'm currently working on the algorithm. Gokhan Dogan discovered a bias in confidence bounds that didn't show up in my original test cases. More examples are helpful. It's possible that the change to the random number generator has already fixed it, but I haven't tested yet. Apparently the old generator wasn't up to the scale of random number usage in MCMC. More examples are always helpful.

Nate Osgood had a paper at the last conference on Bayesian/MCMC methods (I guess we should have written one too, given that we implemented it!).

[gwr]

For completeness:

Lecture on MCMC for System Dynamics by N. Osgood

Pesentation on MCMC for System Dynamics by N. Osgood with examples

2013 Conference Paper on MCMC by N. Osgood

There is also a link to doing sequential Monte Carlo (Particle Filtering) to be found on the YouTube sites above -- albeit using some Russian competitor's software.

And then there is also this great video about sequential monte carlo: The Particel Filter explained without equations

Guido

[WayneZhou2009]

Hi Tom,

I know this is an old post, but I find it very useful to me as someone just beginning to learn how to apply Bayesian approach to calibration of parameters in Vensim. Having carefully read through the above discussion and in particular your replies, I still have some question about creating a policy payoff involving priors.

In your reply above, you gave the code as below ( of course you already clarified that the likelihood should have been log likelihood):

payoff = prior likelihood + point likelihood
prior likelihood = IF THEN ELSE( Time = INITIAL TIME, <expression for parameter likelihood>, 0)
point likelihood = IF THEN ELSE( data var <> :NA:, <expression for data point likelihood - e.g. -((model var - data var)/err sd)^2>, 0)

My question is: for the prior, why is it log likelihood rather than the log pdf? Shouldn't the payoff be like this:

payoff = log ( prior pdf ) + point log likehihood
prior pdf = IF THEN ELSE( Time = INITIAL TIME, <expression for parameter prior pdf>, 0)
point likelihood = same as your original equation

Also, for prior, why Time = Initial time is required?

Many thanks for your time.

[tomfid]

I used likelihood interchangeably with "prior pdf" which is probably a bit more accurate.

The INITIAL TIME check is used because you only compute the prior once, not every time.

The first method of computation, using IF THEN ELSE statements, is obsolete now that you can create a mixed payoff. You probably want:

Code: Select all

*CG
model var|data var/stdDev
*PF
prior likelihood/1

with prior likelihood = log(parameter pdf).

*CG gives you the full likelihood, including a term for the stdDev, so you can estimate that if needed.
*PF or *PI ensures that the prior is computed only once.

[WayneZhou2009]

Hi Tom,

Thank you so much for your reply. In a hypothetical case, such as the rate of identifying some flaws in products, I want to calibrate my model against a set of measurement data whose distribution follows Weibull distribution ("Data Weibull PDF"). So I define my model as a Weibull distribution characterized by a shape parameter and a scale parameter ("Model Weibull PDF"). Then I want to estimate the shape and scale.

Now, suppose prior knowledge and information suggest the Weibull shape parameter should follow a lognormal distribution and the Weibull scale parameter should follow a uniform distribution. So I use lognormal distribution as the prior PDF of shape parameter and uniform distribution as the prior PDF of scale parameter.

In the attached model, I still use the explicit expression of the log likelihood (instead of the mixed payoff that you mentioned), just to help myself better understand the logic behind it.

My payoff (Total log likelihood) is:

-((Model Weibull PDF-Data Weibull PDF)/0.1)^2 + IF THEN ELSE(Time=INITIAL TIME, LN(Prior PDF of scale parameter),0) + IF THEN ELSE(Time=INITIAL TIME, LN(Prior PDF of shape parameter),0)

Please note that the "0.1" is the standard deviation. I just arbitrarily put 0.1 there. Prior PDF of shape parameter is the Lognormal distribution, Prior PDF of scale parameter is the uniform distribution.

My question is that:

(1) If I did this correctly, the result of maximizing the payoff seems to be very sensitive to the standard deviation, which serves as the weight. How to make sure the standard deviation takes a reasonable value?
(2) Is there a need to use the full expression of the likelihood of the Normal distribution, in this case, because of the two LN (prior PDF) in the whole equation?
(3) How to evaluate the role of the prior PDFs (i.e. the 2nd and 3rd components in the additive equation above), given that they might be very weak in some circumstances? In other words, in what circumstances do they carry substantial "value" so that the use of Bayes approach is necessary/meaningful?

Thank you so much for your time!

[tomfid]

A couple things:

In the attached model, I still use the explicit expression of the log likelihood (instead of the mixed payoff that you mentioned), just to help myself better understand the logic behind it.

Often a good idea!

(1) If I did this correctly, the result of maximizing the payoff seems to be very sensitive to the standard deviation, which serves as the weight. How to make sure the standard deviation takes a reasonable value?

I think you're seeing dominance of the scale parameter's contribution via the prior. Setting a larger SD weakens the contribution of the data, so 1/scale in the prior tends to shrink the scale.

(2) Is there a need to use the full expression of the likelihood of the Normal distribution, in this case, because of the two LN (prior PDF) in the whole equation?

No, if you're not estimating the SD. However, you do want to divide the first term (model-data squared) by 2.

(3) How to evaluate the role of the prior PDFs (i.e. the 2nd and 3rd components in the additive equation above), given that they might be very weak in some circumstances? In other words, in what circumstances do they carry substantial "value" so that the use of Bayes approach is necessary/meaningful?

One option would be to vary the SD on the data portion, as you've done. Another would be to add a weight to the priors, so they can be turned down or off as an experiment. It also helps to try some simulations that deviate from the optimum, and watch how the payoff varies - if you have an explicit structure, as here, you can then easily identify whether the curvature is due to the data or the prior.

I think the big issue here is that you're mixing metaphors. The model/data term is evaluating a pdf rather than individual measurements. This means that the measurements are spaced uniformly over the x axis (=time) rather than being distributed with higher frequency at the mode of the pdf. Heuristically, this should still work, but it's more like generalized method of means, and I'm not sure the weighting is valid as a result.

[tomfid]

Try the attached:

Weibull calibration 2.mdl

(4.86 KiB) Downloaded 386 times

The default has 0 weight to priors. I get:

Maximum payoff found at:
Shape parameter = 2.00946
*Scale parameter = 9.82314
Simulations = 70
Pass = 3
Payoff = -289.904

... very close to the true parameters.

If you turn on the prior weight (=1), I get:

Maximum payoff found at:
Shape parameter = 1.98888
*Scale parameter = 9.77772
Simulations = 80
Pass = 3
Payoff = -293.701

Notice that the shape parameter shrinks a bit due to the prior (as you'd expect), but the data still dominates.

Now try a prior with a narrow sigma (=.1) and the wrong value (say, 3 instead of 2):

Maximum payoff found at:
Shape parameter = 2.77034
*Scale parameter = 10.5843
Simulations = 71
Pass = 3
Payoff = -303.959

The best fit is now biased towards the confident prior.

[WayneZhou2009]

Hi Tom, thank you so much for your detailed reply. My further questions are below, highlighted in blue colour:

A couple things:

In the attached model, I still use the explicit expression of the log likelihood (instead of the mixed payoff that you mentioned), just to help myself better understand the logic behind it.

Often a good idea!

(1) If I did this correctly, the result of maximizing the payoff seems to be very sensitive to the standard deviation, which serves as the weight. How to make sure the standard deviation takes a reasonable value?

I think you're seeing dominance of the scale parameter's contribution via the prior. Setting a larger SD weakens the contribution of the data, so 1/scale in the prior tends to shrink the scale.

I used 1/scale as the prior PDF for scale parameter, i.e. a uniform distribution, as a kind of "non-informative prior". By saying "so 1/scale in the prior tends to shrink the scale", did you mean this is a good or bad practice?

(2) Is there a need to use the full expression of the likelihood of the Normal distribution, in this case, because of the two LN (prior PDF) in the whole equation?

No, if you're not estimating the SD. However, you do want to divide the first term (model-data squared) by 2.

First of all, please could you confirm if my understanding is correct: it is assumed that the error (i.e. model - data) follows a normal distribution, and therefore the "standard deviation" in the likelihood function refers to the standard deviation of the error ?

I double-checked the log-likelihood function of normal distribution. The 1st item is a constant which is -n/2*log(2pie), whereas the 2nd item a function of standard deviation, i.e. -n/2*log ((standard deviation)^2). So, if standard deviation is a concern and playes a role in the 2nd item, why should the 2nd item be omitted from the equation?

And, thanks for your reminder about dividing -(model-data)^2 by 2. So, when standard deviation is explicitly considered, this becomes - (model-data)^2/std_dev^2/2

(3) How to evaluate the role of the prior PDFs (i.e. the 2nd and 3rd components in the additive equation above), given that they might be very weak in some circumstances? In other words, in what circumstances do they carry substantial "value" so that the use of Bayes approach is necessary/meaningful?

One option would be to vary the SD on the data portion, as you've done. Another would be to add a weight to the priors, so they can be turned down or off as an experiment. It also helps to try some simulations that deviate from the optimum, and watch how the payoff varies - if you have an explicit structure, as here, you can then easily identify whether the curvature is due to the data or the prior.

I can understand the use of a weight to manipulate the priors. However, for the standard deviation, instead of arbitrarily varying it, would this be a reasonable method - I firstly fit the model without using priors, secondly I can calculate the standard deviation of the errors based on the fitted model, and thirdly I plug the calculated standard deviation into the log likelihood function of the model WITH priors and maximize it. Does this make sense?

I think the big issue here is that you're mixing metaphors. The model/data term is evaluating a pdf rather than individual measurements. This means that the measurements are spaced uniformly over the x axis (=time) rather than being distributed with higher frequency at the mode of the pdf. Heuristically, this should still work, but it's more like generalized method of means, and I'm not sure the weighting is valid as a result.

Thanks for pointing this out. I have realised that my model is a bit confusing and does not serve my original purpose. I will soon upload a new one for further discussion. Thanks.

[WayneZhou2009]

Hi Tom, many thanks for this model and your explanation. It well illustrates the relative role of priors. My question is: what is the consideration behind the formula for Data Weibull, i.e. true scale * ( -LN( RANDOM UNIFORM(0,1,0) ) )^(1/true shape) ?

Try the attached: Weibull calibration 2.mdl

[tomfid]

I think the shrinkage from the non-informative prior is fine.

The standard deviation is indeed the sd of the error.

The log(sd) term can be omitted if you're not estimating the sd in the optimization, because it just shifts the payoff by a constant amount that's invariant with respect to the parameters.

You can iterate to estimate the SD of the errors, but I don't think it makes sense to get the SD of the prior that way. The SD of the prior represents your confidence in the prior, not anything related to the data.

I got the expression for Data Weibull from the wiki page on the distribution as an alternative for generating Weibull vars. The parameterization was a little more straightforward than the Vensim version.

[WayneZhou2009]

Thanks a lot for your quick response Tom. Highly appreciated. The iteration to estimate the SD is for the errors, not for the prior. I did not make myself clear enough. But, the point is, after getting the SD for the errors, it will need to be used for the Total log likelihood which is the sum of the log likelihood and the log priors. And then this Total log likelihood will be maximized in order to get the calibrated shape and scale parameters. This process does not involve any change to the SD of the prior ("sigma" in this case). Did you mean this method does not make sense?

For the log(sd), why is it a constant amount? I think it is dependent upon the parameters. Different combinations of parameters (shape, scale) may lead to different simulation results and thus the SD of error (Model - Data) may be different, and thus the inclusion of log (SD) in the log likelihood payoff may affect the result of maximization to different extents? Sorry if this sounds very trivial and straightforward to you, but I still cannot understand it fully.

Many thanks for your time.

[tomfid]

Iteration to estimate the SD of the errors is fine.

The term log(sd) itself is invariant as long as you aren't changing sd in the optimization. It affects other terms, but that's captured already in other terms.

None of this is ever trivial or easy!