Replicating HGF simulations presented in Mathys et al. 2011 with the current version of the toolbox

[milan-andrejevic]

Hi everyone!

Thanks for the great work making, documenting and improving the toolbox!

In order to gain a better understanding of HGF parameters, I am trying to replicate Figure 10 from the Mathys et al. 2011 Frontiers in Human Neuroscience paper using the current version of the toolbox and ehgf scripts. I am having some trouble getting the third level mu parameter (log volatility of a tendency) to rise as it does in Figure 10, so I thought I'd ask for advice!

Here I will make a reproducible example with notes. First I generate a sequence of binary outcomes (0 and 1) similar to the one in Figure 10.

% The sequence is composed of three segments:
% 1) 100 pseudo randomly generated trials for which the probability of occurrence of 1 is set to .9, and the probability of occurrence of 0 is set to .1
SeqL = randsrc(100, 1, [0,1; .1, .9]);

% 2) 9 blocks of 20 trials each, alternating between the .9 probability for 0, and .9 probability for 1. 
SeqS = [randsrc(20,1, [0,1; .9, .1]); randsrc(20,1, [0,1; .1, .9]);...
        randsrc(20,1, [0,1; .9, .1]); randsrc(20,1, [0,1; .1, .9]);...
        randsrc(20,1, [0,1; .9, .1]); randsrc(20,1, [0,1; .1, .9]);...
        randsrc(20,1, [0,1; .9, .1]); randsrc(20,1, [0,1; .1, .9]);...
        randsrc(20,1, [0,1; .9, .1])];
% 3) a copy of sequence generated in 1)
u=[SeqL; SeqS; SeqL];

Then I simulate the model relying on information from the figure notes (please do let me if I am making any errors here). The notes say: theta = .2, omega2 = -2.3, kappa2 = 1.6 Instead of theta, which is no longer in the config scripts in the new versions of HGF, we are specifying omega3. Since theta=exp(omega3) we can calculate omega3 as omega3=ln(theta)=-1.6094.

In the figure we can see that starting values for mu (mu0) were set to 0 for the first level, -2 for the second level, and 0 for the third level. I also assume that sigma2 and 2 were set to 0.

sim = tapas_simModel(u,...
                     'tapas_ehgf_binary',...
                      [NaN -2 0, NaN 0 0, NaN 0 0, 1 1.6, NaN -2.3 log(.2)],... 
                      ...%mu1 2 3,sgm1 2 3,rh1 2 3, k1 2, om0, om2, om3]
                     'tapas_unitsq_sgm');

tapas_hgf_binary_plotTraj(sim)

Here's what we get:

As you can see the mu3 shifts around but doesn't generally increase/accumulate after a volatile period of alternating blocks. I thought my omega3 parameter was off, perhaps too high so I thoroughly tested shifting it around.

In fact I've tried playing with each parameter, shifting one around while I keep others to the default I outlined above, and also tried using the regular 'tapas_hgf_binary' script instead of ehgf, but I couldn't get closer to reproducing the figure, which makes me think I must be making an error somewhere or there may be something more I need to tweak?

Any advice would be highly appreciated!!

Best wishes, Milan

[chmathys]

Hi Milan,

You can reproduce Figure 10 when using the original tapas_hgf_binary (not the new tapas_ehgf_binary). The input stimulus sequence from the paper is here: u.mat.zip

Code:

u = load(u.mat)

sim = tapas_simModel(u,...
                     'tapas_hgf_binary',...
                     [NaN -1.8 0 NaN 1 1 NaN 0 0 1 1.6 NaN -2.3 -1.61],...
                     'tapas_unitsq_sgm',...
                     5,...
                     123456789);

tapas_hgf_binary_plotTraj(sim)

This pretty much exactly reproduces the figure:

Note that the figure in the paper was made about two years before there was an HGF toolbox.

[chmathys]

Having said this, I'd also remark that if you're interested in input sequences that drive the volatility estimate steeply upwards, you should use inputs on a continuous scale if you can.

It's quite hard for any filter to infer changes in volatility on top of changes in contingency from binary inputs alone.