Kruschke chapter 10

[Dirk Nachbar]

I was trying to implement the DBDA p278 model, below is my attempt, stan doesn't like the m parameter (as it's int). How would I get around this? 

Are there any other examples of model probability a la Kruschke?

data {

 int N;

 int y[N];

}

parameters {

 simplex[2] mPriorProb; 

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

}

transformed parameters {

 real kappa;

 int<lower=1, upper=2> m[N]; #????

 simplex[2] omega;

 kappa <- 12;

 omega[1] <- .25;

 omega[1] <- .75;

}

model {

 y ~ bernoulli(theta);

 theta ~ beta(omega[m]*(kappa-2)+1, (1-omega[m])*(kappa-2)+1);

 m ~ categorical(mPriorProb);

}

Dirk

[Mike Lawrence]

I believe you need to put the initialization of the variable m in the parameters section, not the transformed parameters section. 

Also, kappa & omega look like data variables to me, as they are always the same fixed value you're assigning into them.

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

--
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+unsubscribe@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Bob Carpenter]

First, in Stan, transformed parameters must not only be declared (as you
have done), but also defined (given values, which you didn't do) in the
transformed parameter block.

You are correct that Stan does not support integer parameters.
You have to marginalize them out of the density (you can reconstruct
them if you want, but you lose information from sampling relative
to just working in expectation).

In your example, the m parameter provides mixture "responsibilities".
There's a chapter in the manual on mixture models that explains how to
implement them by marginalizing out the discrete parameter. There's
also a chapter on latent discrete parameters that has more general
examples.

Given that you define all your transformed parameters directly, they
should really be defined in the transformed data block.

If you can't figure it out after reading the above and giving it
a try, let us know and we can help further.

What we should really do is just translate all the examples in his
book into Stan to make this easy for people.

- 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.

[Bob Carpenter]

> On Aug 11, 2016, at 4:00 PM, Mike Lawrence <mike....@gmail.com> wrote:
>
> I believe you need to put the initialization of the variable m in the parameters section, not the transformed parameters section.

No. Stan won't let you set anything in the parameters block, nor
will it let you define integer parameters.

> Also, kappa & omega look like data variables to me, as they are always the same fixed value you're assigning into them.

Absolutely. It's more efficient that way, though in this trivial
model, it won't matter.

- Bob

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>

> --
> 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.

[Mike Lawrence]

> I believe you need to put the initialization of the variable m in the parameters section, not the transformed parameters section.

No.  Stan won't let you set anything in the parameters block, nor
will it let you define integer parameters.

Ah, by "initialize" I meant declared (i.e. "int<lower=1, upper=2> m[N];"), but wasn't aware of the no-ints limitation (though makes sense given my vague understanding of the transforms that happen behind the scenes).

[Mike Lawrence]

The 2nd edition of that book has stan code (https://sites.google.com/site/doingbayesiandataanalysis/software-installation), though possibly not for every example?

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--
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+unsubscribe@googlegroups.com.

[Bob Carpenter]

Two replies in one.

> The 2nd edition of that book has stan code
> (https://sites.google.com/site/doingbayesiandataanalysis/software-installation),
> though possibly not for every example?

There's a chapter on Stan, I don't think every example
got translated.

> Ah, by "initialize" I meant declared (i.e. "int<lower=1, upper=2> m[N];"),

We also allow initialization of parameters through the
interfaces. The declaration is just that --- a declaration
of the type. It doesn't do anything like initialization to
a concrete value.

> but wasn't aware of the no-ints limitation (though makes sense given my vague understanding of the transforms that happen behind the scenes).

HMC is the restriction---it only handles continuous
parameters. Systems like PyMC3 that do HMC and discrete
parameters alternate discrete sampling and HMC.

- Bob

[Avraham Adler]

On Thursday, August 11, 2016 at 9:46:41 AM UTC-4, Dirk Nachbar wrote:

I was trying to implement the DBDA p278 model, below is my attempt, stan doesn't like the m parameter (as it's int). How would I get around this? 

[snip]
 

Dirk

Kruschke, on page 415 of the second edition (subtitled A tutorial with R, JAGS, and Stan) when discussing Stan's limitation of not having discrete parameters, specifically gives the example of chapter 10.3.2 (the one on page 278) as something that cannot directly be translated from JAGS to Stan (second paragraph of section 14.5).

Avi

[Bob Carpenter]

I'll translate this model to Stan if you can
send me the JAGS code. Or tell me the file name
in the DBDA download---I don't have the book.

Basically any time you see a simplex[2], you know
this isn't a good coding because it can be replaced
with a real<lower=0, upper=1> (the first component
of the simplex) for much more straightforward models.

- Bob


> On Aug 11, 2016, at 3:46 PM, 'Dirk Nachbar' via Stan users mailing list <stan-...@googlegroups.com> wrote:
>

> --
> 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.

[Dirk Nachbar]

Thanks Bob, here is the Jags

model {

for ( i in 1:N ) {

y[i] ˜ dbern( theta )

}

theta ˜ dbeta( omega[m]*(kappa-2)+1 , (1-omega[m])*(kappa-2)+1 )

omega[1] <- .25

omega[2] <- .75

kappa <- 12

m ˜ dcat( mPriorProb[] )

mPriorProb[1] <- .5

mPriorProb[2] <- .5

[Andrew Gelman]

It would be good to clear this up.  Once we have it running, we can post it on blog and send an email to Kruschke so he can fix this in the next edition of his book!

A

[Avraham Adler]

He doesn't say it's impossible, but says it may require "clever programming contrivances" :)

On Thursday, August 11, 2016 at 4:55:39 PM UTC-4, Andrew Gelman wrote:

It would be good to clear this up.  Once we have it running, we can post it on blog and send an email to Kruschke so he can fix this in the next edition of his book!

A

On Aug 11, 2016, at 10:55 AM, Avraham Adler <avraha...@gmail.com> wrote:

On Thursday, August 11, 2016 at 9:46:41 AM UTC-4, Dirk Nachbar wrote:

I was trying to implement the DBDA p278 model, below is my attempt, stan doesn't like the m parameter (as it's int). How would I get around this? 

[snip]
 

Dirk

Kruschke, on page 415 of the second edition (subtitled A tutorial with R, JAGS, and Stan) when discussing Stan's limitation of not having discrete parameters, specifically gives the example of chapter 10.3.2 (the one on page 278) as something that cannot directly be translated from JAGS to Stan (second paragraph of section 14.5).

Avi

--
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+unsubscribe@googlegroups.com.

[Andrew Gelman]

That's just wrong.  From my perspective, the mixture model as written in Stan is more direct and less of a "contrivance" than the latent-variable expression in Jags.

A

[Bob Carpenter]

The original JAGS code from Kruschke is in

archive: DBDA2Eprograms.zip
file: Jags-Ydich-Xnom1subj-MbernBetaModelComp.R
URL: https://sites.google.com/site/doingbayesiandataanalysis/software-installation

Here's the JAGS model:

model {
for ( i in 1:N ) {

y[i] ~ dbern( theta )
}
theta ~ dbeta( omega[m]*(kappa-2)+1 , (1-omega[m])*(kappa-2)+1 )
omega[1] <- .25

omega[2] <- .75
kappa <- 12

m ~ dcat( mPriorProb[] )

mPriorProb[1] <- .5
mPriorProb[2] <- .5
}

And here's the Stan version (which just evaluates all those
constant expressions and plugs them in):

data {
int<lower=0> N;
int<lower=0, upper=1> y[N];
}
parameters {
real<lower=0, upper=1> theta;
}
model {
y ~ bernoulli(theta);
target += log_mix(0.5, beta_lpdf(theta | 2.75, 8.25),
beta_lpdf(theta | 8.25, 2.75));
}

As seen from Gelman's and Kruschke's comments, "contrivance" is often
in the eye of the beholder. I'm with Kruschke on this one (though I
wouldn't have written the JAGS code with all those constants) but
that may be because I learned Bayesian stats through BUGS.

Here's the fit from RStan 2.11.1 (remember, don't use 2.10!)
using all default settings (the first three lines are from Krushcke's
JAGS example in the file linked above):

> N=9
> z=6
> y = c(rep(0, N - z), rep(1, z));
> fit <- stan("kruschke-db.stan", data=c("y", "N"));
> print(fit)

Inference for Stan model: kruschke-db.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
theta 0.67 0.00 0.14 0.32 0.6 0.70 0.77 0.87 1409 1
lp__ -7.69 0.02 0.77 -9.69 -8.0 -7.38 -7.14 -7.09 2089 1

For comparison, here's the fit you get with a uniform prior on
theta, also using all default settings:

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
theta 0.64 0.00 0.14 0.34 0.55 0.64 0.74 0.88 2639 1

- Bob

P.S. Here are the Stan files if you want to run yourself:

[Christian Schuhegger]

Sorry for reopening such an old topic, but I have a question around your example. The Kruschke chapter is about model comparison and the valuable information is in the 'm' parameter, which basically tells you, which model to prefer.

How would I recover the 'm' information from the Stan model you've given?