preventing underflow in a mixture model

Martin Šmíra

unread,

Oct 2, 2014, 6:55:33 PM10/2/14

to stan-...@googlegroups.com

Hello Stan,

I have a question about a variation of a reinforcement learning model (for Iowa gambling task), which I already discussed here: https://groups.google.com/forum/?hl=cs#!topic/stan-users/TjY3wQqUS2g

I want to extend the model with a mixture of a learning rule. First, I created an individual model, where the mixing determined for each subject s seperately. So choice[s,t] for each subject is determined either by Delta learning rule or Decay learning rule for all trials t.

This model actually works as expected (see below).

transformed parameters {

for (s in 1:n_s) {

for (t in 1:(n_t - 1)) {

lp_tmp_1[t] <- categorical_logit_log(choice[s,t+1], Ev_delta * theta);

lp_tmp_2[t] <- categorical_logit_log(choice[s,t+1], Ev_decay * theta);

}

lp_parts[s,1] <- log(mix) + sum(lp_tmp_1);

lp_parts[s,2] <- log1m(mix) + sum(lp_tmp_2);

}

model {

for (s in 1:n_s)

increment_log_prob(log_sum_exp(lp_parts[s]));

}

Problem starts when I change assumtion of the model. Now I want the mixture of the learning rule to be determined for whole group of subjects at once. So that all subjects choices are following either Delta or Decay rule. The model is below.

transformed parameters {

for (s in 1:n_s) {

for (t in 1:(n_t - 1)) {

lp_tmp[1,s,t] <- categorical_logit_log(choice[s,t+1], Ev_delta * theta);

lp_tmp[2,s,t] <- categorical_logit_log(choice[s,t+1], Ev_decay * theta);

}

lp_parts[1] <- log(mix) + sum(lp_tmp[1]);

lp_parts[2] <- log1m(mix) + sum(lp_tmp[2]);

}

model {

for (s in 1:n_s)

increment_log_prob(log_sum_exp(lp_parts));

}

generated quantities {

prob <- softmax(lp_parts - max_lp);

z <- bernoulli_rng(prob[1]);

}

Finally getting to the problem. There is 100 trials for 22 subjects. Therefore accumulated log probability in variable lp_parts is getting very low (reaching minus thousands). So my humble guess is an underflow. Practically when I run the model with more than a few subjects then it results in no mixing. Chains for parameter z (mixing indicator) are stuck on either 0.00 or 1.00 after a few jumps on warmup phase. To be clear, in a one run several chains are on 0 and several on 1.

I've tried to reformulate the computation into:

model {

increment_log_prob(log_sum_exp(lp_parts) + max_lp);

}

, where max_lp is max(lp_parts).

OR:

model {

increment_log_prob(log1p(exp(min_lp - max_lp)) + max_lp);

}

, where max_lp is max(lp_parts) and min_lp is min(lp_parts)

But it didn't help much.

Full model is to be found here: https://github.com/neekro/Stan_translate/blob/master/MODEL-MIX_pvl_vs_pvl_d.stan

And full print from stanfit is here: https://github.com/neekro/Stan_translate/blob/master/model-mix_output.txt

You can see that min_lp and max_lp can still differ by almost 2000.

Is there a way?

Cheers,

Martin

Martin Šmíra

unread,

Oct 11, 2014, 8:41:14 PM10/11/14

to stan-...@googlegroups.com

I couldn't find a way overcome this issue yet.

Still hoping for an answer.

Thanks.

Bob Carpenter

unread,

Oct 12, 2014, 1:39:09 PM10/12/14

to stan-...@googlegroups.com

Thanks for the reminder --- I accidentally marked the mail as
read and never came back to it.

I'll take a look and get back to you as soon as I can.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

Bob Carpenter

unread,

Oct 14, 2014, 12:14:41 PM10/14/14

to stan-...@googlegroups.com

I'm afraid I still don't understand the model well enough
from the code fragments.

Is there a writeup somewhere of the model you're trying to fit?

Also, you can print with more digits to see if results
are truly pinned.

You won't be able to make log_sum_exp() any more arithmetically
stable. It's already implemented so that

m <- max(a,b);
log_sum_exp(a,b) = m + log(exp(a - m) + exp(b - m))

Problems with gradients can also cause issues with not being able to move.

If you use softmax explicitly, that can easily underflow in
some dimensions because it puts the result on the linear scale.
So it's more stable to use categorical_logit.

- Bob

Yes, if the log prob is minus thousands, then the linear prob
is zero. But it still shouldn't underflow if you're using
categorical_logit or softmax.

> Practically when I run the model with more than a few subjects then it results in no mixing. Chains for parameter z (mixing indicator) are stuck on either 0.00 or 1.00 after a few jumps on warmup phase. To be clear, in a one run several chains are on 0 and several on 1.
>
> I've tried to reformulate the computation into:
>
> model {
> increment_log_prob(log_sum_exp(lp_parts) + max_lp);
> }
> , where max_lp is max(lp_parts).
>
> OR:
>
> model {
> increment_log_prob(log1p(exp(min_lp - max_lp)) + max_lp);
> }
> , where max_lp is max(lp_parts) and min_lp is min(lp_parts)
>
> But it didn't help much.
>
> Full model is to be found here: https://github.com/neekro/Stan_translate/blob/master/MODEL-MIX_pvl_vs_pvl_d.stan
> And full print from stanfit is here: https://github.com/neekro/Stan_translate/blob/master/model-mix_output.txt
>
> You can see that min_lp and max_lp can still differ by almost 2000.
>
> Is there a way?
>
> Cheers,
> Martin
>
>
>

Martin Šmíra

unread,

Oct 22, 2014, 8:56:30 AM10/22/14

to stan-...@googlegroups.com

Hi Bob,

In this article: Steingroever, H., Wetzels, R., & Wagenmakers, E.-J. (2013). Validating the PVL-Delta model. Frontiers in Decision Neuroscience, 4: 898.

(https://dl.dropboxusercontent.com/u/12798592/Steingroever_Frontiers_2013.pdf)

In this article they are using identical model, but wihout the extenstion I'm interested in. In section 2.2 the model (PVL-Delta) is explained and later in section 4.1.1 is figure with the graphical model and its description.

All I want to do is to assess whether Delta learning rule used in PVL-Delta model fits the data better then alternative model PVL (without Delta) which uses Decay learning rule. And I want to assess it by setting the model as a mixture model where the learning rules will be mixed. From posterior mixing proportion of the model I can compute BF indicating wich learning rule (or model) is more appropriate for my data. This procedure is called Product space method.

(Lodewyckx, T., Kim, W., Lee, M. D., Tuerlinckx, F., Kuppens, P., & Wagenmakers, E.-J. (2011). A tutorial on Bayes factor estimation with the product space method. Journal of Mathematical Psychology, 55(5), 331–347. doi:10.1016/j.jmp.2011.06.001)

I think that my implementation of the model is correct, but as number of subjects gets higher (even 5 subjects can be too much) then the model performs bad mixing. Log probabilities of each mixture part is reaching minus thousands on the log scale. When I compute posterior mixing probability (variable prob) and look at its traceplot then I can see that the values are jumping between extremes and spend just a little time in the midlle values (eg. 0.5) and its histogram is therefore bimodal. This happens in the case where I have just about 5 subjects in the group, if I have more then there is no mixing at all (some chains on 0 and other on 1).

Here is my model: https://github.com/neekro/Stan_translate/blob/master/MODEL-MIX_pvl_vs_pvl_d.stan

How do I actually find out if I have a problem with gradient? I tried the grad_test, but I don't understand the output much. I read on user list that there shouldn't be unsensible values, but it is not the case here. Only thing that surprised me was that Log probability for each chain has varied from -695472 to -3085.15. See below full report of one chain.

TEST GRADIENT MODE

Log probability=-29109.7

param idx value model finite diff error

0 -0.075878 -3936.73 -3936.73 3.36836e-006

1 -1.38366 -1601.93 -1601.93 8.24962e-006

2 1.59553 -2393.25 -2393.25 5.27943e-006

3 -0.0187773 -64300.1 -64300.1 6.00796e-006

4 0.959427 -950.028 -950.028 1.13291e-007

5 -1.28679 -344.207 -344.207 5.10236e-006

6 1.23755 557.387 557.387 -5.28166e-006

7 -1.93995 -7319.34 -7319.34 1.62422e-006

8 0.836111 -2709.14 -2709.14 1.62192e-006

9 0.774026 -68.9066 -68.9066 -1.12684e-006

10 1.23084 -93.7432 -93.7432 6.65539e-007

11 -0.75303 -3.0343 -3.0343 1.54631e-006

12 0.00912743 -0.0877146 -0.0877153 6.7612e-007

13 1.85108 4.57733 4.57733 4.02398e-007

14 -1.53884 10.3859 10.3859 6.34157e-007

15 1.34729 -54.8447 -54.8447 -9.04885e-007

16 -1.0827 -27.9604 -27.9604 1.84978e-006

17 -0.203475 31.0896 31.0896 -2.86708e-006

18 -1.60571 50.5896 50.5896 -3.43172e-006

19 -1.63125 7.13502 7.13501 2.5133e-006

20 1.97514 -16.4945 -16.4945 2.7873e-006

21 -0.670087 52.5602 52.5602 5.91858e-007

22 0.888881 -421.857 -421.857 -2.47141e-006

23 0.0664499 15.6815 15.6815 6.22744e-007

24 1.62147 -55.9499 -55.9499 1.47571e-006

25 0.12502 -413.708 -413.708 -9.61516e-007

26 0.467332 -593.541 -593.541 -2.63238e-006

27 1.95469 2.2516 2.2516 8.36731e-009

28 -1.57313 3.27532 3.27533 -2.80828e-007

29 -1.93522 10.0557 10.0557 2.13032e-006

30 0.393893 -272.975 -272.975 1.18992e-006

31 0.245291 -19.8614 -19.8614 4.50591e-007

32 -1.61851 -324.621 -324.621 -3.13852e-006

33 -1.46912 -73.3448 -73.3448 -8.28408e-007

34 1.2659 -2.57249 -2.57249 -8.8145e-007

35 1.87396 -2.00954 -2.00954 -5.63978e-007

36 -0.416298 -14.3928 -14.3928 1.60834e-006

37 -1.87633 -269.758 -269.758 -6.11611e-007

38 -1.06442 -590.657 -590.657 1.47104e-006

39 -1.47521 -185.624 -185.624 2.29841e-006

40 -0.257608 -39.2774 -39.2774 -2.19497e-006

41 -0.874811 -198.09 -198.09 6.46929e-008

42 1.43053 -2.53264 -2.53264 -1.70653e-006

43 -1.39745 1.48499 1.48499 -3.36794e-007

44 -0.805542 -322.605 -322.605 -3.44079e-007

45 1.51495 -1.62916 -1.62916 -7.59306e-007

46 0.253802 -77.0643 -77.0643 -5.27784e-006

47 0.957385 -12.0946 -12.0946 4.15508e-008

48 0.270758 -98.753 -98.753 -3.66487e-007

49 1.23026 -3.68422 -3.68422 2.07055e-007

50 -0.57321 -267.153 -267.153 2.99208e-006

51 0.687278 -1.87706 -1.87706 1.42161e-006

52 -0.78214 -73.2005 -73.2005 6.86488e-007

53 1.81516 -146.137 -146.137 -2.33733e-006

54 -0.327891 -9.68709 -9.68709 6.47486e-007

55 -1.73519 3.18269 3.18268 1.85553e-006

56 0.732026 -22.3775 -22.3775 4.41305e-007

57 0.0647289 -25.7771 -25.7771 -8.72306e-009

58 -1.95623 -0.477764 -0.477759 -4.46529e-006

59 -1.10183 -3.88192 -3.88192 2.09369e-006

60 0.498988 -13.3387 -13.3387 1.43029e-006

61 0.386984 -11.43 -11.43 -7.05476e-007

62 -0.479387 14.0488 14.0488 -8.28766e-007

63 -1.63466 -4.71691 -4.71691 -1.06316e-007

64 1.88064 -71.7772 -71.7772 -1.52674e-006

65 -1.404 2.71395 2.71395 7.03046e-007

66 0.539788 -23.1226 -23.1226 -1.88553e-006

67 0.0352478 -17.6737 -17.6737 2.58643e-006

68 -0.471258 -20.4139 -20.4139 3.92704e-006

69 0.259487 -11.2613 -11.2613 4.18764e-007

70 0.938896 -77.34 -77.34 2.25564e-007

71 0.862999 -24.4834 -24.4834 -3.57391e-006

72 0.325911 14.4106 14.4106 1.94015e-006

73 0.991226 2.87349 2.87349 1.25122e-007

74 1.9448 -2169.03 -2169.03 -1.66934e-006

75 1.52298 -654.536 -654.536 1.16155e-006

76 -0.324174 -395.181 -395.181 -8.28778e-007

77 -0.639936 -270.971 -270.971 8.71003e-007

78 0.267065 -410.198 -410.198 2.76365e-006

79 1.18812 -553.672 -553.672 1.12984e-006

80 -0.509854 -204.734 -204.734 -2.36759e-006

81 -0.0825078 -295.477 -295.477 -1.49933e-006

82 1.92485 -824.309 -824.309 -1.72193e-006

83 -0.0287695 -305.686 -305.686 2.17355e-006

84 0.593136 -617.85 -617.85 1.21059e-006

85 0.0382916 -424.06 -424.06 -1.15963e-006

86 -0.379679 -403.907 -403.907 -3.04161e-006

87 -1.71943 -173.106 -173.106 2.06056e-007

88 -0.50263 -545.465 -545.465 4.85116e-007

89 -1.15066 -166.113 -166.113 -6.74951e-007

90 0.483965 -707.953 -707.953 -1.45403e-006

91 -1.04656 -537.104 -537.104 -6.28818e-007

92 0.548405 -824.352 -824.352 -5.32883e-006

93 1.75991 -470.169 -470.169 1.3625e-006

94 1.11494 -966.023 -966.023 -3.46528e-007

95 -1.65596 -202.765 -202.765 -1.03903e-006

Cheers,

Martin

Bob Carpenter

unread,

Oct 22, 2014, 4:33:45 PM10/22/14

to stan-...@googlegroups.com

The gradient looks fine. The new Stan 2.5 manual has an explanation of
what's reported. It's basically a comparison of the built-in autodiff
and finite differences on the unconstrained scale.

We don't really do Bayes factor calculations because Andrew prefers
model expansion, so I haven't had any experience doing this.

That product-space paper costs $36 until I get back inside
Columbia's network. I found something similar in terms of authors
and content online on Michael Lee's site:

http://www.socsci.uci.edu/~mdlee/LodewyckxEtAl2009.pdf

It looks like from figure 3 that the mixture is at the individual level,
but share a mixture parameter hat-pi (which confused me, because hats are
usually used for estimates, not naming parameters).

Do both models fit the data the way you'd expect without
the mixture?

Can you compute Bayes factors for other models using the same
technique?

I'm not very familiar with these calculations, but I'd watch out
for normalization issues --- the densities in Stan are not
normalized in the sampling formulation, only in the increment_log_prob
with the _log density functions.

- Bob

Martin Šmíra

unread,

Oct 23, 2014, 7:56:19 PM10/23/14

to stan-...@googlegroups.com

Hi Bob,

great, I see that Stan 2.5 is out now, I'm happy about the progress.

I thought that the paper is freely available, can you access it here (below)?

https://ppw.kuleuven.be/okp/_pdf/Lodewyckx2011ATOBF.pdf

The Bayes Factor is not really important here. Only thing where the model differs from the original one is the possibility of sampling from categorical_logit with two different parameters. First is givven by Delta learning rule and another by Decay leraning rule:

// Delta Learning Rule

Ev_delta[choice[s,t]] <- (1 - a_ind[s]) * Ev_delta[choice[s,t]] + a_ind[s] * v;

// Decay Lerning Rule

Ev_decay <- Ev_decay * a_ind[s];

Ev_decay[choice[s,t]] <- Ev_decay[choice[s,t]] + v;

choice[s,t+1] ~ categorical_logit(Ev_delta * theta);

choice[s,t+1] ~ categorical_logit(Ev_decay * theta);

If this would mix properly, then I can easily get the BF from the posterior mixing proportion.

Q; Do both models fit the data the way you'd expect without
the mixture?
- Yes, I did, but I better check it again.

Q: Can you compute Bayes factors for other models using the same technique?

- Yes, I've already used this method, actually on a similar (PVL-Delta) model. This model was comparing two groups of subjects whether they come from the same hierarchical group distribution. And it worked like a charm, actually just in Stan, Jags version performed miserably (poor mixing). Actually I've heard that Jags (BUGS also I guess) in general has problem with this method. I guess that this approach is computationally demanding. Especially in this case where probabilities of each mixture part is on such a small scale.

What do you mean by the normalization issues? I use "~" sampling statements only in 8 cases and they all are just normal(0, 1).

Thanks a lot for looking at it!

Martin

Bob Carpenter

unread,

Oct 25, 2014, 1:23:44 PM10/25/14

to stan-...@googlegroups.com

On Oct 23, 2014, at 7:56 PM, Martin Šmíra <smira....@gmail.com> wrote:

> The Bayes Factor is not really important here. Only thing where the model differs from the original one is the possibility of sampling from categorical_logit with two different parameters. First is givven by Delta learning rule and another by Decay leraning rule:
>
> // Delta Learning Rule
> Ev_delta[choice[s,t]] <- (1 - a_ind[s]) * Ev_delta[choice[s,t]] + a_ind[s] * v;
> // Decay Lerning Rule
> Ev_decay <- Ev_decay * a_ind[s];
> Ev_decay[choice[s,t]] <- Ev_decay[choice[s,t]] + v;
>
> choice[s,t+1] ~ categorical_logit(Ev_delta * theta);
> {{OR (as a mixture)}}
> choice[s,t+1] ~ categorical_logit(Ev_decay * theta);

I'm confused by pretty much all of this. If you want
to code a mixture, you want to use

increment_log_prob(log_sum_exp(lp1, lp2));

where lp1 and lp2 are the log probs of the two components being mixed,
including the mixing proportions. For example, a mixture that I'd write
in math as:

p(y) = lambda * normal(y|mu1,sigma1) + (1 - lambda) * normal(y|mu2,sigma2)

is coded as:

increment_log_prob(log_sum_exp(log(lambda) + normal_log(y,mu1,sigma1),
log1m(lambda) + normal_log(y,mu2,sigma2)));

First, it looks like a_ind[s] is some kind of mixing proportion. But
it looks like you're mixing some kind of scalar multiplier? I don't see
why the mixture defining Ev_delta is done on the logit scale. But maybe a_ind
isn't a mixing proportion?

And then there's a move from t to t+1 that I don't get.

Why would you be multiplying by theta in both categorical logits? What
is the type of theta and why is it the same in both cases? And where does
it make sense to multiply on the logit scale? That has the effect of
flattening or attenuating the probabilities depending on whether the
multipler is below or above zero.

If you could write down in ordinary math notation what the distribution
is supposed to be for choice[s,t] (I also don't get the t vs. t+1 above),
then it should be easy to code in Stan.

> If this would mix properly, then I can easily get the BF from the posterior mixing proportion.
>
> Q; Do both models fit the data the way you'd expect without
> the mixture?
> - Yes, I did, but I better check it again.
>
> Q: Can you compute Bayes factors for other models using the same technique?
> - Yes, I've already used this method, actually on a similar (PVL-Delta) model. This model was comparing two groups of subjects whether they come from the same hierarchical group distribution. And it worked like a charm, actually just in Stan, Jags version performed miserably (poor mixing). Actually I've heard that Jags (BUGS also I guess) in general has problem with this method. I guess that this approach is computationally demanding. Especially in this case where probabilities of each mixture part is on such a small scale.

If you have the A/B models in Stan and JAGS, could you send them along?
We're looking to write some papers on comparisons and looking for examples
that behave differently in Gibbs vs. HMC/NUTS.

> What do you mean by the normalization issues? I use "~" sampling statements only in 8 cases and they all are just normal(0, 1).

Categorical logit is what's being mixed, and it will always be normalized,
so that's not an issue.

> I thought that the paper is freely available, can you access it here (below)?
> https://ppw.kuleuven.be/okp/_pdf/Lodewyckx2011ATOBF.pdf

The link is to an illegally pirated copy that someone's hosting
on a web page. As much as I hate closed-source journals, I'm not
advocating stealing. The journal does let you buy open access for your
paper (their site says fees range from US$500 to US$5000:

http://www.elsevier.com/about/open-access/open-access-journals

but only two authors seem to have taken them up on it in the last six
years:

http://www.journals.elsevier.com/journal-of-mathematical-psychology/open-access-articles/

- Bob

Martin Šmíra

unread,

Oct 25, 2014, 5:52:53 PM10/25/14

to stan-...@googlegroups.com

Hello Bob,

I've posted a link on github to the full model two posts ago. In the previous post I just wanted to make it clearer about the mixture, but nevermind. I'm posting the full model below:

# PVL-Delta in Stan

data {

int<lower=1> n_s; // Total # subjects

int<lower=1> n_t; // # trials

int<lower=0,upper=4> choice[n_s, n_t]; // # subj. x # trials matrix with

// choices

real<lower=-35,upper=2> net[n_s, n_t]; // Net amount of wins + losses

// (# subj. x # trials matrix)

real<lower=0,upper=1> mix; // Prior mixing proportion

}

parameters {

// Group-level mean parameters

real mu_A_pr;

real mu_a_pr;

real mu_w_pr;

real mu_c_pr;

// Group-level standard deviation

real<lower=0,upper=1.5> sd_A;

real<lower=0,upper=1.5> sd_a;

real<lower=0,upper=1.5> sd_w;

real<lower=0,upper=1.5> sd_c;

vector[n_s] raw[4];

}

transformed parameters {

vector[2] lp_parts;

matrix[n_s,n_t-1] lp_tmp[2];

// Individual-level paramters

vector<lower=0,upper=1>[n_s] A_ind;

vector<lower=0,upper=5>[n_s] w_ind;

vector<lower=0,upper=1>[n_s] a_ind;

vector<lower=0,upper=5>[n_s] c_ind;

for (s in 1:n_s) {

A_ind[s] <- Phi(mu_A_pr + sd_A * raw[1,s]);

w_ind[s] <- Phi(mu_w_pr + sd_w * raw[2,s]) * 5;

a_ind[s] <- Phi(mu_a_pr + sd_a * raw[3,s]);

c_ind[s] <- Phi(mu_c_pr + sd_c * raw[4,s]) * 5;

}

for (s in 1:n_s) { // loop over subjects

real theta;

vector[4] Ev_delta;

vector[4] Ev_decay;

// Trial 1

Ev_delta <- rep_vector(0, 4); // assign 0 to all elements of Ev_delta

Ev_decay <- rep_vector(0, 4); // assign 0 to all elements of Ev_delta

theta <- 3 ^ c_ind[s] - 1;

// Remaining trials

for (t in 1:(n_t - 1)) {

real v;

if (net[s,t] >= 0)

v <- (net[s,t]) ^ A_ind[s];

else

v <- -w_ind[s] * fabs(net[s,t]) ^ A_ind[s];

// Delta Learning Rule

Ev_delta[choice[s,t]] <- (1 - a_ind[s]) * Ev_delta[choice[s,t]] + a_ind[s] * v;

// Decay Lerning Rule

Ev_decay <- Ev_decay * a_ind[s];

Ev_decay[choice[s,t]] <- Ev_decay[choice[s,t]] + v;

lp_tmp[1,s,t] <- categorical_logit_log(choice[s,t+1], Ev_delta * theta);

lp_tmp[2,s,t] <- categorical_logit_log(choice[s,t+1], Ev_decay * theta);

}

lp_parts[1] <- log(mix) + sum(lp_tmp[1]);

lp_parts[2] <- log1m(mix) + sum(lp_tmp[2]);

}

model {

# Prior on the group-level mean parameters

# probit scale [-Inf, Inf]

mu_A_pr ~ normal(0, 1);

mu_a_pr ~ normal(0, 1);

mu_w_pr ~ normal(0, 1);

mu_c_pr ~ normal(0, 1);

for (i in 1:4)

raw[i] ~ normal(0, 1);

increment_log_prob(log_sum_exp(lp_parts));

}

generated quantities {

real<lower=0,upper=1> mu_A;

real<lower=0,upper=5> mu_w;

real<lower=0,upper=1> mu_a;

real<lower=0,upper=5> mu_c;

vector<lower=0,upper=1>[2] prob;

int z; // mixing variable

mu_A <- Phi(mu_A_pr);

mu_w <- Phi(mu_w_pr) * 5;

mu_a <- Phi(mu_a_pr);

mu_c <- Phi(mu_c_pr) * 5;

prob <- softmax(lp_parts);

z <- bernoulli_rng(prob[1]);

}

I have two formulas on how to compute "expectancy valence" (Ev). So I created two variables Ev_delta and Ev_decay. You can see their definition in the mode, Ev_delta is just on one line, but Ev_decay is on two lines - one line is updating all (4) elements in the vector, and second one affects just one element.

I want to know which formula for computing Ev is better fitting my data. And thats what I want to achieve with the mixture model.

I used categorical_logit() here as an equivalent to categorical(softmax()) function, is that correct? So originaly it looked like this:

prob <- softmax(Ev_delta * theta);

choice[s,t+1] <- categorical(prob);

Where theta [0, 5] is a free parameter used for either flattening or attenuating the probabilities (as you said). And that has some desired implications for the model. And it is a vector of same size as number of subjects.

To be more specific, The model is used for a task (Iowa Gabling Task). There you have 4 decks of cards and some decks are more rewarding then others. In each trial a subject will choose a card and receive some monetary reward or penalty. In the model variable Ev_delta (or Ev_decay) is a vector of size 4 and it's purpose is to grasp subject's expected gains from each deck. At the beginning all elements of Ev are assigned 0 and after each trial, depending on reward, the Ev is updated by either Delta rule or by Decay rule. In practise only one these rules is used (I want both!). After that this Ev is transformed to probabilities (with a little help of theta) and used as a prediction of a next deck choice (choice[s,t+s]).

"t+1" is used there because we want to skip the first choice[s,1], which is just random.

Anyway, I just run the models without the mixture part and it all works well (both for Ev_delta and Ev_decay). We've used this model (without mixture) many times so I'm quite sure that the problem comes with the mixture implementation.

Do you still want me to write down the distribution for choice[s,t] in math notation? I'm not sure if I would work it out.

I'll send you those models you wanted for comparison later. I don't have the Jags one on my disk, but I'll try to get it.

Cheers,

Martin

Herra Huu

unread,

Oct 25, 2014, 9:43:18 PM10/25/14

to

Open access article about the model: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3786232/

By simply looking at the model definitions, I can see at least one possible source of problems. Why parameter a_ind is used in both decay and delta learning rules in the following code?

// Delta Learning Rule

Ev_delta[choice[s,t]] <- (1 - a_ind[s]) * Ev_delta[choice[s,t]]

+ a_ind[s] * v;

// Decay Lerning Rule

Ev_decay <- Ev_decay * a_ind[s];

Ev_decay[choice[s,t]] <- Ev_decay[choice[s,t]] + v;

Shouldn't there be separate parameters for both rules?

Martin Šmíra

unread,

Oct 26, 2014, 2:53:11 PM10/26/14

to stan-...@googlegroups.com

Thank you for your suggestion Herra,

but I've already tried a model version with two separate a_ind parameters some time ago (and now again). It didn't help. Some chains are stuck on one mixture part, other chains are stuck on the other.

But I was thinking that it shouldn't matter if there is one a_ind parameter or two in this model. Since there are only two options in the model - either all subjects are using Delta rule (first mixture) or the Decay rule (second mixture). And therefore it shouldn't lead to any interaction between a_ind information from those two mixture parts. Right?

I've also tried to remove group-level hierarchical distributions for a_ind, w_ind, A_ind and c_ind. This also didn't help.

Bob Carpenter

unread,

Oct 26, 2014, 3:32:54 PM10/26/14

to stan-...@googlegroups.com

I think there's two things you want to do to get a top-level
mixture here.

First, you need to separately parameterize each model.
That means, for instance, breaking theta into two variables,
say theta_delta and theta_decay. Similarly, you'll need
v_delta and v_decay because that looks like a shared parameter.
If you have the individual models, it should be easy to figure
out what needs to be duplicated. It looks like there's quite
a few.

Second, I think you want to declare mix as real<lower=0,upper=1>,
and then estimate it. You can then assign bernoulli_rng(mix) to
z in each iteration; but it's more efficient to do inference directly
with mix rather than draws z.

- Bob

Bob Carpenter

unread,

Oct 26, 2014, 3:41:45 PM10/26/14

to stan-...@googlegroups.com

Oops -- didn't see Herra's comments before responding.

Yes, it does matter. You need to split the parameters.
You can share the data. Both sub-models are evaluated at each step
and then mixed.

The density that you are then defining is:

p(lambda, theta1, theta2, y)

= lambda * M1(theta1, y) + (1 - lambda) * M2(theta2, y)

or on the log scale:

log p(lambda, theta1, theta2, y)

= log_sum_exp(log(lambda) + log M1(theta1, y),
log1m(lambda) + log M2(theta2, y))

where theta1 and theta2 are disjoint parameters and y is data,
which does not have to be duplicated.

I just wanted to make sure this was the density you intended
to estimate. If so, if you split the parameters and estimate
lambda (what you called mix), then you should be good to go.

- Bob

Herra Huu

unread,

Oct 26, 2014, 6:25:46 PM10/26/14

to

Shouldn't the whole joint density p(m, theta_1, theta_2, y) be the following:

p(m=1) * p(theta_1|m=1) * p(y|theta_1,m=1) * p(theta_2|m=1) +
p(m=2) * p(theta_2|m=2) * p(y|theta_2,m=2) * p(theta_1|m=2)
, where p(theta_1|m=2) and p(theta_2|m=1) are arbitrarily chosen pseudo priors/linking densities.

But the density in the current code, as far as I can see, is more something like this (not quite, because of the a_ind etc):
[p(m=1) * p(y|theta_1,m=1) + p(m=2) * p(y|theta_2,m=2)] *
p(theta_2|m=2) * p(theta_1|m=1)

Is there any examples of Stan implementation of product-space method?

Martin Šmíra

unread,

Oct 27, 2014, 5:38:27 PM10/27/14

to stan-...@googlegroups.com

I've tried your suggestions and I also removed group-level parameter to simplify the model (see below). This model is not working properly either.

# PVL-Delta in Stan

data {

int<lower=1> n_s; // Total # subjects

int<lower=1> n_t; // # trials

int<lower=0,upper=4> choice[n_s, n_t]; // # subj. x # trials matrix with

// choices

real<lower=-35,upper=2> net[n_s, n_t]; // Net amount of wins + losses

// (# subj. x # trials matrix)

}

parameters {

// Individual-level paramters

vector<lower=0,upper=1>[n_s] A_ind_delta;

vector<lower=0,upper=5>[n_s] w_ind_delta;

vector<lower=0,upper=1>[n_s] a_ind_delta;

vector<lower=0,upper=5>[n_s] c_ind_delta;

vector<lower=0,upper=1>[n_s] A_ind_decay;

vector<lower=0,upper=5>[n_s] w_ind_decay;

vector<lower=0,upper=1>[n_s] a_ind_decay;

vector<lower=0,upper=5>[n_s] c_ind_decay;

real<lower=0,upper=1> mix; // mixing proportion

}

model {

vector[2] lp_parts;

matrix[n_s,n_t-1] lp_tmp_delta;

matrix[n_s,n_t-1] lp_tmp_decay;

for (s in 1:n_s) { // loop over subjects

real theta_delta;

real theta_decay;

vector[4] Ev_delta;

vector[4] Ev_decay;

// Trial 1

Ev_delta <- rep_vector(0, 4); // assign 0 to all elements of Ev_delta

Ev_decay <- rep_vector(0, 4); // assign 0 to all elements of Ev_delta

theta_delta <- 3 ^ c_ind_delta[s] - 1;

theta_decay <- 3 ^ c_ind_decay[s] - 1;

// Remaining trials

for (t in 1:(n_t - 1)) {

real v_delta;

real v_decay;

if (net[s,t] >= 0) {

v_delta <- net[s,t] ^ A_ind_delta[s];

v_decay <- net[s,t] ^ A_ind_decay[s];

} else {

v_delta <- w_ind_delta[s] * fabs(net[s,t]) ^ A_ind_delta[s];

v_decay <- w_ind_decay[s] * fabs(net[s,t]) ^ A_ind_decay[s];

}

// Delta Learning Rule

Ev_delta[choice[s,t]] <- (1 - a_ind_delta[s]) * Ev_delta[choice[s,t]]

+ a_ind_delta[s] * v_delta;

// Decay Lerning Rule

Ev_decay <- Ev_decay * a_ind_decay[s];

Ev_decay[choice[s,t]] <- Ev_decay[choice[s,t]] + v_decay;

lp_tmp_delta[s,t] <- categorical_logit_log(choice[s,t+1], Ev_delta * theta_delta);

lp_tmp_decay[s,t] <- categorical_logit_log(choice[s,t+1], Ev_decay * theta_decay);

}

increment_log_prob(log_sum_exp(log(mix) + sum(lp_tmp_delta),

log1m(mix) + sum(lp_tmp_decay)));

}

Is it really possible to estimate the mix parameter this way? I thought that the mix parameter won't be updated in this case.

I've generated synthetic data, one batch generated with Delta learning rule and another batch generated with Delta learning rule. And I've run the model each data separately. Unfortunately, the model didn't manage to identify the data correctly. The estimation of parameter mix was on every run around 0.35 with deviation around 0.25. It seems that the mix parameter is not updated properly. The model always seems to prefer Decay rule. In a different model implementation I observed that the log probability for model 1 (Delta) was reaching - 100 000. In this run I estimated just 6 subjects, part of the output is below:

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

mix 0.350 0.015 0.212 0.015 0.182 0.338 0.492 0.833 200 1.004

lp_parts[1] -108710.229 11570.920 140596.043 -540786.316 -119906.849 -61629.303 -26066.002 -7624.266 148 0.998

lp_parts[2] -196.001 0.892 4.376 -204.451 -198.515 -195.709 -193.026 -189.887 24 1.071

lp__ -271.403 1.196 6.709 -284.466 -274.895 -270.913 -266.875 -261.058 31 1.102

I've checked the model at least 10 times find my mistake, but I couldn't.

Herra: I've already implemented product space method in Stan. I hope that correctly. Here is the code: https://github.com/neekro/Stan_translate/blob/master/pvl_d_mix_one_1a_(1-std).stan

It is also the PVL-Delta model, but the product space method was used to infer whether two groups of subjects differ on A parameter.

Thanks,

Martin

Bob Carpenter

unread,

Oct 27, 2014, 6:07:44 PM10/27/14

to stan-...@googlegroups.com

I've tried to understand the tutorial on Bayes factor product
calculations.

1. Lodewyckx, Lee and Wagenmakers.
A general computational method for estimating Bayes factors with the
product space method.
preprint version: http://www.socsci.uci.edu/~mdlee/LodewyckxEtAl2009.pdf
The official version is in the Journal of Mathematical Psychology from 2011
and a bunch more authors.

2. BUGS examples, vol2, Pines
http://www.medicine.mcgill.ca/epidemiology/joseph/courses/EPIB-682/pines.pdf

There's also a paper by Geyer with "umbrella sampling" in the title in
the MCMC handbook, but it doesn't add much to the above discussions.

I think there may be trouble with the way the factor method
approximates reversible jump MCMC by using Gibbs and pseudo-priors.
I'm not sure it'll work in HMC. I haven't been able to get it to work,
anyway. And it looks awfully fiddly even if it could be made to work.

If someone could get the BUGS vol2/pines example working, that'd point
the way to doing these models in general.

I'll bring the topic up at the Stan meeting tomorrow if there's time,
but we already have a busy agenda this week.

- Bob

Reply all

Reply to author

Forward