problem with a simple correlation model

[Martin Šmíra]

I'm again using this book and trying to implement their Bugs models to Stan: https://webfiles.uci.edu/mdlee/BB_Free.pdf. The model is on page 63. Graphical model and Bugs model code can be found there. 

I think that I don't have any mistake in the model. But the sampling gets too autocorrelated and the samples also seems to be a bit different on the average. Look at the attached pdf. Negative correlation is more probable in Stan compared to Jags (same with WinBugs) and the difference is also visible in higher values of "r". 

I have tried many things to solve it, but the results were always unstable and if stable, then a bit different from Jags.

It starts to be problematic when I input higher "sigmaerror". When I put = c(.03,1), sampling is neat.

What I've noticed is that sigma[2] is twice higher on average in Stan and is having problem reaching lower values. 

Any idea what goes wrong? I never understand why these autocorrelations are happening and what triggers it. 

Stan Code:  

# clears workspace: 

rm(list=ls(all=TRUE)) 

library(rstan)

model <- "

# Pearson Correlation With Uncertainty in Measurement

data{ 

  int<lower=0> n;

  vector[2] x[n];

  vector[2] sigmaerror;

}

parameters{

  vector[2] mu;

  vector<lower=0>[2] lambda;

  real<lower=-1,upper=1> r;

  vector[2] y[n];

} 

transformed parameters{

  vector<lower=0>[2] sigma;

  cov_matrix[2] T;

  #Reparameterization

  sigma[1] <- 1 / sqrt(lambda[1]);

  sigma[2] <- 1 / sqrt(lambda[2]);

  T[1,1] <- 1 / lambda[1];

  T[1,2] <- r * sigma[1] * sigma[2];

  T[2,1] <- r * sigma[1] * sigma[2];

  T[2,2] <- 1 / lambda[2];

}

model{

  # Priors

  mu ~ normal(0,sqrt(1000));

  lambda ~ gamma(.001,.001);

  r ~ uniform(-1,1);

  # Data

  for (i in 1:n){

    y[i] ~ multi_normal(mu,T);

    for (j in 1:2){

      x[i,j] ~ normal(y[i,j],sigmaerror[j]);

    }

  }

}"

x <- matrix(c(.8, 102,

              1,  98,

             .5, 100,

             .9, 105,

             .7, 103, 

             .4, 110,

            1.2,  99,

            1.4,  87,

             .6, 113, 

            1.1,  89,

            1.3,  93), nrow=11,ncol=2,byrow=T) 

n <- nrow(x) # number of people/units measured

sigmaerror = c(0.03, 10)

data <- list(x=x, n=n, sigmaerror=sigmaerror) # to be passed on to Stan

samples <- stan(model_code = model, data = data,# init=myinits,

               iter = 5000, chains = 3)           

print (samples)

Thank you for any help.

Martin

[Andrew Gelman]

Martin:

I'm not sure, but these dgamma(.001,.001) priors can create serious problems; see here:

http://www.stat.columbia.edu/~gelman/research/published/taumain.pdf

I don't know what's going on in this particular example, but this is not a model I'd ever want to be fitting.

Do you get the same problem if you use uniform priors on sigma1 and sigma2?  Also, for this model, we might prefer the LKJ parameterization of the covariance matrix.  Stan has various covariance models so you can model that matrix directly, no need for the T[1,1], T[2,1], etc. stuff.

Andrew

P.S.  I'm also ccing EJ because I'm sure he would want his examples to work in Stan also.

--
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/groups/opt_out.
<compare.pdf>

[Martin Šmíra]

Yes, I've tried using uniform prior for sigma1 (0,2) and sigma2 (0,50), but it doesn't help. 

Do you think LKJ parameterization would help? I'll try it tomorrow, seems complicated. 

I'm actually doing it for EJ. I'm working on implementing all those models from the book. 

Martin

[Bob Carpenter]

On 9/7/13 3:44 PM, Andrew Gelman wrote:
> ...

> Also, for this model, we might prefer the
> LKJ parameterization of the covariance matrix. Stan has various covariance models so you can model that matrix
> directly, no need for the T[1,1], T[2,1], etc. stuff.

We might want to break out the small size cases in
those definitions to make them more efficient.

Martin -- the LKJ models are easy. The best thing to
do is use a correlation matrix type and prior, then
scale it the same way you're doing for the 2 x 2 case.

- Bob

[Ben Goodrich]

The LKJ stuff is easy. But in the 2x2 case, an LKJ prior on the correlation matrix with eta = 1 is the same as a uniform prior on (-1,1) for the correlation. So, in this case, I think its all how you model the precisions / variances / standard deviations.

Ben

[Bob Carpenter]

On 9/7/13 7:00 PM, Martin �m�ra wrote:
> Yes, I've tried using uniform prior for sigma1 (0,2) and sigma2 (0,50), but it doesn't help.
> Do you think LKJ parameterization would help? I'll try it tomorrow, seems complicated.

I don't think it'll help control the sampling. But you'll
want to do something like this for higher dimensions than 2.

I reparameterized your model with a uniform prior for scale
instead of your vague inverse gamma for precision.

I also rewrote it to be a bit more efficient. It may not be
ideal for pedagogical purposes --- we struggle with the balance
between writing clear and writing fast models (ideally, Stan would be
smart enough to turn the former into the latter!).

It looks like the model mixes well and produces posterior draws for
r (correlation) that look like the values I see in figure 5.4.
And it's recovering reasonable mu values. I'm not seeing the
wild multimodal posterior densities you plotted.

Where did your inputs come from and did they match the inputs
that generated fig 5.4 in the notes? From the text,
it looks like

sigmaerror = c(0.03,1)

was the input. But in your code, it was

sigmaerror = c(0.03,10).

Even if I use more error, I recover roughly the same values, though effective
sample size is quite a bit lower per iteration.

You also have this:

mu ~ normal(0,sqrt(1000));

which is an informative prior on the scale you're at,
because the posterior mean of mu[2] is 98, whereas sqrt(1000)
is only 32, and there are only 12 noisy data points!

For instance, changing it to

mu ~ normal(0,sqrt(10000));

raises the posterior mean estimate of mu[2] from 98 to 100.

This prior may also not be interacting well with the noise model that has
scale sd = 10 in the more extreme parameterization.

In general, it can be problematic when you're estimating values
near boundaries. What happens in Stan is that we logit transform
bounded values (and then scale and translate them). Values near
boundaries on (-1,1) represent logit values near -infinity or +infinity,
which can be difficult numerically.

Also, with the two different scales, Stan has to do more work.
Ideally, the data itself would be a bit more standardized.

> I'm actually doing it for EJ. I'm working on implementing all those models from the book.

:-) Small world, of course.

- Bob

P.S. I see you have kappa coefficients coming up next. I'd suggest
looking at something like Dawid and Skene's 1979 measurement error
model for data annotation in addition to kappa.

[Martin Šmíra]

Thank you all for advices.  

Bob:

> Where did your inputs come from and did they match the inputs 
> that generated fig 5.4 in the notes?  From the text, 

> it looks like  sigmaerror = c(0.03,1) .

 

You're right, the figure 5.4 is with sigmaerror <- .03,1. With these sigmaerror values everything goes well. But as soon as you change it to sigmaerror <- .03,10, as is required in exercise 5.2.2 (in the book), then it gets more interesting. See the pictures below, where I compared posteriors from Jags, Bugs and our 2 Stan models. First pair is around 10k samples with several chains. The second one is about a half of previous.

About your code Bob. What is " ' " supposed to do there y[i]' ~ multi_normal(mu,T);? And is it better to use col(x,j) instead of having array of vectors? I thought that in manual the second option was prefered.

I've also tried mu ~ normal(0,sqrt(10000)), but it's resulting in similar correlation posterior. 

Can I actually build on results from Jags and BUGS? Or they can be wrong in this case?

Yep, seems like Baysian world is quite small :-)

Martin

[Eric-Jan Wagenmakers]

Dear Stan users,

Thanks for the feedback. First off, yes, when you assume large measurement error the convergence for sigma becomes problematic; in fact, we have started to look in to this just last week (a coincidence; Martin was unaware of this). The convergence problems may arise because of model misspecification (the data were not generated with large measurement error; instead, the were generated with no measurement error, and then the amount of measurement error was simply assumed to be present). The convergence problem is also present in WinBUGS.

Thanks for all the other suggestions for improvement

Cheers,
E.J.

[Sergio Polini]

Il 08/09/2013 01:00, Martin �m�ra ha scritto:
> I'm actually doing it for EJ. I'm working on implementing all those
> models from the book.

I'm a learner, I'm trying to rewrite from BUGS to STAN the models from
the ARM book, so I've looked at the book by Lee & Wagenmakers too.

The example about Cha Sa-soon (�5.5) is really amusing, but it could
also be interesting (and eventually not so easy to rewrite) because:
a) one can't pass a vector such as z <- c(rep(NA, 949), 30) to Stan;
b) the censored data can be treated as missing data (Stan Manual, �9.2,
p. 44), but one can't define parameters such as

int<lower=15,upper=25> z[N_cens];

because:
(i) the error message I got said that integer parameters are not allowed;
(ii) the binomial() function requires an integer lhs.

So after looking at "Integrating out Censored Values" (Manual, p. 45)
I've tried:

data{
int<lower=0> N_cens; // 949
int<lower=0> n; // 50
int<lower=0> z_obs; // 30
}
parameters {
real<lower=0.25,upper=1.0> p;
}
model {
z_obs ~ binomial(n, p);
lp__ <- lp__ + N_cens *
log(binomial_cdf(25, n, p) - binomial_cdf(14, n, p));
}

My results look reassuring:

> print(chasasoon.sf, probs=c(0.025,0.5,0.975), digits=2)
...
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
p 0.40 0.00 0.0 0.39 0.40 0.41 1871 1
lp__ -152.26 0.03 0.7 -154.62 -151.97 -151.75 489 1

I hope this helps... if I'm right.
I hope someone helps me if I'm wrong ;-)

Sergio

[Bob Carpenter]

On 9/8/13 4:26 PM, Sergio Polini wrote:

> Il 08/09/2013 01:00, Martin Šmíra ha scritto:
>> I'm actually doing it for EJ. I'm working on implementing all those
>> models from the book.
>
> I'm a learner, I'm trying to rewrite from BUGS to STAN the models from the ARM book, so I've looked at the book by Lee &
> Wagenmakers too.
>

> The example about Cha Sa-soon (§5.5) is really amusing, but it could also be interesting (and eventually not so easy to

> rewrite) because:
> a) one can't pass a vector such as z <- c(rep(NA, 949), 30) to Stan;

Right. No NA in Stan. Our data types are based
on C++ double and int.

> b) the censored data can be treated as missing data (Stan Manual, §9.2, p. 44), but one can't define parameters such as

>
> int<lower=15,upper=25> z[N_cens];
>
> because:
> (i) the error message I got said that integer parameters are not allowed;

> (ii) the binomial() function requires an integer lhs.

Right.

> So after looking at "Integrating out Censored Values" (Manual, p. 45) I've tried:
>
> data{
> int<lower=0> N_cens; // 949
> int<lower=0> n; // 50
> int<lower=0> z_obs; // 30
> }
> parameters {
> real<lower=0.25,upper=1.0> p;
> }
> model {
> z_obs ~ binomial(n, p);
> lp__ <- lp__ + N_cens *
> log(binomial_cdf(25, n, p) - binomial_cdf(14, n, p));
> }

Cool --- that's one observation of 30 and N_cens observations of an
unknown value in a known range.

We really need to get the following functions for each density
FOO built in:

g_FOO(x,y,theta) = log(FOO_cdf(x,theta) - FOO_cdf(y,theta))

We use it in actual truncation models, too, and it's not
the most arithmetically stable way to implement the function,
nor the most efficient given the shared parameters theta and
the use of auto-diff to do all the derivatives.

> My results look reassuring:
>
> > print(chasasoon.sf, probs=c(0.025,0.5,0.975), digits=2)
> ...
> mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
> p 0.40 0.00 0.0 0.39 0.40 0.41 1871 1
> lp__ -152.26 0.03 0.7 -154.62 -151.97 -151.75 489 1
>

Nice.

> I hope this helps... if I'm right.
> I hope someone helps me if I'm wrong ;-)

Looks good to me.

- Bob

[Bob Carpenter]

On 9/8/13 8:49 AM, Martin �m�ra wrote:
> Thank you all for advices.
>
> Bob:
>
> > Where did your inputs come from and did they match the inputs
> > that generated fig 5.4 in the notes? From the text,
> > it looks like sigmaerror = c(0.03,1) .
>
> You're right, the figure 5.4 is with sigmaerror <- .03,1. With these sigmaerror values everything goes well. But as soon
> as you change it to sigmaerror <- .03,10, as is required in exercise 5.2.2 (in the book), then it gets more interesting.
> See the pictures below, where I compared posteriors from Jags, Bugs and our 2 Stan models. First pair is around 10k
> samples with several chains. The second one is about a half of previous.

I didn't have a problem fitting my model with sigmaerror <- 0.03,10.

> <https://lh6.googleusercontent.com/-MFVr9AaiGDE/UixpCiCPMDI/AAAAAAAAAKM/QuhNWFosi-k/s1600/big_compare2.jpeg>
>
>
> <https://lh6.googleusercontent.com/-HDFKxGLUuEg/Uixo8eIA1NI/AAAAAAAAAKE/N2DP5_1X5xE/s1600/big_compare.jpeg>

>
> About your code Bob. What is " ' " supposed to do there y[i]' ~ multi_normal(mu,T);?

Vector transpose -- it turns a row vector, which is
what y[i] picks out, to a column vector (or just "vector"),
which is what multi_normal expects.

> And is it better to use col(x,j)
> instead of having array of vectors? I thought that in manual the second option was prefered.

So yes, if you made y an array of vectors, that'd be faster than
using y[i]'.

> I've also tried mu ~ normal(0,sqrt(10000)), but it's resulting in similar correlation posterior.
>
> Can I actually build on results from Jags and BUGS? Or they can be wrong in this case?

Unfortunately, all these methods can fall over in the same situations,
which typically involve highly variable scales in the posterior. What happens
is that there are sharp transitions in log posterior (for instance, when a hierarchical
variance parameter goes to zero), which are hard for all of these methods to
deal with. RMHMC will be better, but it can be slow per iteration and
limited in scalability.

I'm worried about the two versions of this model being different in Stan. If
everything else agreed other than my model, my model might be wrong. Seems
more likely than all the other ones being wrong.

- Bob

[Martin Šmíra]

Bob,

the posterior of your model was actually different just because of your prior on sigma(unif(0,20)). I made some new plots :)  

Both pairs differs just with prior on sigma, either uniform(0,20) or invGamma(.001,.001). Both trials are using 4 chains, each 10k iterations - 2k for warmup. I specified initials only for WinBUGS.

Seems like Stan doesn't handle gamma distribution well. And your model seems to fit better. 

Can I ignore the warning about matrix being positive definite? Or is it better to give bounds to r (-0.9999,0.9999)?

Thanks for your time Bob.

Martin

[Bob Carpenter]

On 9/9/13 8:44 AM, Martin �m�ra wrote:
> Bob,
>
> the posterior of your model was actually different just because of your prior on sigma(unif(0,20)). I made some new
> plots :)
> Both pairs differs just with prior on sigma, either uniform(0,20) or invGamma(.001,.001). Both trials are using 4
> chains, each 10k iterations - 2k for warmup. I specified initials only for WinBUGS.
> Seems like Stan doesn't handle gamma distribution well. And your model seems to fit better.
>

> <https://lh4.googleusercontent.com/-CP1L271Q-HY/Ui27vIJiG-I/AAAAAAAAAKc/lUigBO8KKFo/s1600/gamma_unif1.jpeg>
>
>
> <https://lh3.googleusercontent.com/-kO-G1OKX6RE/Ui27xUaIY1I/AAAAAAAAAKk/vBACGmDPZww/s1600/gamma_unif2.jpeg>

>
>
> Can I ignore the warning about matrix being positive definite?

Yes, the problem arises during the simulation of the
position of the parameters under Hamiltonian dynamics.
We need to trade off the time granularity of the simulation
(step size) vs. number of steps. With smaller step size,
we get more accurate transition approximations, but it
takes more steps. This balance is tuned for optimal
sampling during warmup. The resulting balance occassionally
rejects because of too-large steps using only gradient
taking you outside of support for constrained variables
like covariance matrices.

> Or is it better to give bounds to r (-0.9999,0.9999)?

That will probably do it, but whether it's better or not is a matter
of where you think the results will be. If you get estimates near the
boundary, then you want wider bounds (closer to 0 and 1). The other
alternative is to use a stronger prior pushign values to 0.
That's why we like the LKJ priors --- you can set up a correlation matrix
that favors the identity with an adjustable strength, then use a scale
vector (i.e., the sd of the dimensions) to convert it to a correlation
matrix. So you independently control amount of shrinkage toward
identity correlation and the scales of the dimensions.

> Thanks for your time Bob.

You're welcome. The examples you've brought up have been fun and
we want to make sure Stan is usable in as many situations and for
as many kinds of models as possible.

Speaking of which, do you have plans to distribute the models when
you're done coding them? Or do a comparison with JAGS/BUGS? That'd
be fantastic, I think, for us and the BUGS and JAGS developers.
What we're interested in is for efficiency and scalability are:

(a) speed to convergence, especially from diffuse
initializations,

(b) mixing speed after convergence in terms of number
of effective samples per second, and

(c) memory use.

We're interested not only in averages, but also variation, though I'm
not so sure how to measure that other than histograms of the measurements
above. We've been finding that our chains can be pretty variable in how long they take.

- Bob

[Martin Šmíra]

Yes , the book is full of fun and interesting examples. I recommend to look further into it.

All the examples in the book are already modeled in WinBUGS and Jags and are prepared to run from both R and Matlab.

Stan will be hopefully the next option for the readers. There is not a plan to do a direct comparison, but I think that the code will be avaliable for this purpouse. 

But I'm free to hire to work on that! :)

Martin

[Martin Šmíra]

Hi Sergio,

 I didn't have this one implemented to Stan, I'll use it, thank you for that.

But!  Why did you put "14" in " - binomial_cdf(14, n, p)"? There should be "15". 

Btw: Sampling for this model is like 100x faster compared to Jags nad Bugs.

Martin

[Martin Šmíra]

Bob C. :

 

P.S.  I see you have kappa coefficients coming up next.  I'd suggest
       looking at something like Dawid and Skene's 1979 measurement error
       model for data annotation in addition to kappa.

 You mentioned this kappa example, I already implemented this one to Stan. I did it in a simple way using the same code as is in the EJ's book and I didn't have any problem with that model. 

Martin

[Bob Carpenter]

I just meant in the order of the book. It's great
that there's a reference to Bayesian kappa coefficient
estimates --- I've been suggesting to people that they
do that for years.

I work on data annotation/rating/coding tasks, so
evaluating rater agreement is a problem near and dear
to my heart. There are some really nifty Bayesian
models in the area (Girolami and crew have one for ordinals
with Gaussian process priors; Raykar et al. have one that
includes an automatically trained system as a rater.

It has some psychological relevance because it's all
about humans making categorical judgments and their
accuracies and biases --- so it's also related to education
models like IRT if you start thinking about things like
item difficulty.

- Bob

[Sergio Polini]

Il 11/09/2013 03:11, Martin �m�ra ha scritto:
> Hi Sergio,
>
> I didn't have this one implemented to Stan, I'll use it, thank you for
> that.

It was an (amusing, and) interesting model that helped me to acquaint
myself with "lp__ tricks".

> But! Why did you put "14" in "- binomial_cdf(14, n, p)"? There should
> be "15".

I suppose that binomial_cdf() in Stan is similar to pbinom() in R, i.e.

binomial_cdf(x, size, p) = pbinon(x, size, p) = P[X <= x]

and P[15 <= X <= 25] = P[X <= 25] - P[X < 15] = P[X <= 25] - P[X <= 14].
In R:

> p <- 0
> for (i in 15:25) p <- p + choose(50,i) * 0.4^i * (1-0.4)^(50-i)
> p
[1] 0.8887012
> pbinom(25, 50, 0.4) - pbinom(14, 50, 0.4)
[1] 0.8887012

Sergio

[Martin Šmíra]

Sergio, 

this makes sense, but I think that Stan actually compute it in a different way than R. I think that Stan does some kind of discrete transformation. Is it right Bob? 

Try run your model without the observed trial. You'll get samples of p which are slightly below 0.4. That don't make sense. Above that, I was comparing posterior of theta (p) from Bugs, Jags and Stan (in my version with 15). See below. Ignore BUGS (green) density, there seems to be some issue. 

But in logic. If we have censored observations scoring from 15 to 25. Expected value will be 20 (out of 50) -> expected theta (p) is 0.4. With one more observed trial, which is scoring above expectation should mode theta a bit higher. Right?

[sp_r...@yahoo.it]

Il giorno mercoledì 11 settembre 2013 15:34:15 UTC+2, Martin Šmíra ha scritto:

Sergio, 

this makes sense, but I think that Stan actually compute it in a different way than R. I think that Stan does some kind of discrete transformation. Is it right Bob? 

I hope that Stan *_cdf() functions are what they are said to be, i.e. cumulative distribution functions ;-)

Try run your model without the observed trial. You'll get samples of p which are slightly below 0.4. That don't make sense. Above that, I was comparing posterior of theta (p) from Bugs, Jags and Stan (in my version with 15). See below. Ignore BUGS (green) density, there seems to be some issue. 

I think that exact numbers like "0.4" do not exist in a probabilistic world, even less in a sampled world ;-)
However, if I run the Bugs model from the book I get:

> mean(theta)
[1] 0.3987868
> median(theta)
[1] 0.3988205

These numbers are slightly below 0.4 too, but I'd say that they are the reference numbers.
When I run my model I get:

> mean(p_full)
[1] 0.3989067
> median(p_full)
[1] 0.3989678

and I can't see any substantial differences.

If I run my model without the observed trial, I get:

> mean(p_949)
[1] 0.3981429
> median(p_949)
[1] 0.3981435

"a bit" lower...
 

But in logic. If we have censored observations scoring from 15 to 25. Expected value will be 20 (out of 50) -> expected theta (p) is 0.4. With one more observed trial, which is scoring above expectation should mode theta a bit higher. Right?

How much is "a bit"?
I've tried a crude measure:

> bit <- mean(c(rep(0.4, 949),0.6)) - 0.4
> bit
[1] 0.0002105263

That is, I think that you should mistrust a (sampled) distribution whose mean or median is about 0.41.

IMHO, of course.

Sergio

[Bob Carpenter]

On 9/11/13 1:13 PM, sp_r...@yahoo.it wrote:
> Il giorno mercoledì 11 settembre 2013 15:34:15 UTC+2, Martin Šmíra ha scritto:
>

> Sergio,
> this makes sense, but I think that Stan actually compute it in a different way than R. I think that Stan does some
> kind of discrete transformation. Is it right Bob?
>
>

> I hope that Stan *_cdf() functions are what they are said to be, i.e. cumulative distribution functions ;-)

Yes, they are. And the truncation just calls them in
the obvious way (described in the manual).

> Try run your model without the observed trial. You'll get samples of p which are slightly below 0.4. That don't make
> sense. Above that, I was comparing posterior of theta (p) from Bugs, Jags and Stan (in my version with 15). See
> below. Ignore BUGS (green) density, there seems to be some issue.
>
>

> I think that exact numbers like "0.4" do not exist in a probabilistic world, even less in a sampled world ;-)

Right --- any point has zero density in a continuous distribution.

> However, if I run the Bugs model from the book I get:
> |
> > mean(theta)
> [1] 0.3987868
> > median(theta)
> [1] 0.3988205
> |
>
> These numbers are slightly below 0.4 too, but I'd say that they are the reference numbers.
> When I run my model I get:
> |
> > mean(p_full)
> [1] 0.3989067
> > median(p_full)
> [1] 0.3989678
> |
>
> and I can't see any substantial differences.
>
> If I run my model without the observed trial, I get:
> |
> > mean(p_949)
> [1] 0.3981429
> > median(p_949)
> [1] 0.3981435
> |
>
> "a bit" lower...

Is that surprising?

>
> But in logic. If we have censored observations scoring from 15 to 25. Expected value will be 20 (out of 50) ->

Only if the mean is 20 --- that is if your censoring is centered.
It's usually not in most realistic cases.

> expected theta (p) is 0.4. With one more observed trial, which is scoring above expectation should mode theta a bit
> higher. Right?

I would think so. What's the MCMC error at this point?
Do you get the same result with different inits?

> How much is "a bit"?
> I've tried a crude measure:
> |

> > bit <- mean(c(rep(0.4,949),0.6)) - 0.4
> > bit
> [1] 0.0002105263
> |
>
> That is, I think that you should mistrust a (sampled) distribution whose mean or median is about 4.1.
>
> IMHO, of course.

Never trust anything!

- Bob

[Bob Carpenter]

On 9/11/13 9:34 AM, Martin �m�ra wrote:
> Sergio,
> this makes sense, but I think that Stan actually compute it in a different way than R. I think that Stan does some kind
> of discrete transformation. Is it right Bob?

We don't have any discrete transformations, but I'm
not sure what you're talking about transforming here.

For positive constrained parameters, we log transform them,
but then adjust for the Jacobian of the inverse transfrom so
we get the right sampling behavior.

> Try run your model without the observed trial. You'll get samples of p which are slightly below 0.4. That don't make
> sense. Above that, I was comparing posterior of theta (p) from Bugs, Jags and Stan (in my version with 15). See below.
> Ignore BUGS (green) density, there seems to be some issue.
>
> But in logic. If we have censored observations scoring from 15 to 25. Expected value will be 20 (out of 50) -> expected
> theta (p) is 0.4. With one more observed trial, which is scoring above expectation should mode theta a bit higher. Right?
>

> <https://lh5.googleusercontent.com/-EoOcH4g6RFo/UjBvjupWAoI/AAAAAAAAAK0/7dxnu_vs4xQ/s1600/compare_using_Sergios_model2.jpeg>

If those are the posteriors, it looks like there's something
wrong with the BUGS model --- or at least JAGS and Stan
agree.

- Bob

[Bob Carpenter]

On 9/11/13 2:26 AM, Sergio Polini wrote:

> Il 11/09/2013 03:11, Martin Šmíra ha scritto:
>> Hi Sergio,
>>
>> I didn't have this one implemented to Stan, I'll use it, thank you for
>> that.
>
> It was an (amusing, and) interesting model that helped me to acquaint myself with "lp__ tricks".
>
>> But! Why did you put "14" in "- binomial_cdf(14, n, p)"? There should
>> be "15".
>
> I suppose that binomial_cdf() in Stan is similar to pbinom() in R, i.e.

Yes, it's the usual the cumulative distribution function
for the distribution.

http://en.wikipedia.org/wiki/Cumulative_distribution_function

We need to write that out more clearly in the manual. I'll
put it on the to-do list. The functional forms need to go in, too,
but they're very hairy, so we haven't gotten around to it yet.

We're adding cdf() and ccdf() in Stan 2 on the log scales (ccdf() is
1 - cdf()). And we'll be putting in ongoing work to make them more efficient
and more arithmetically stable.

> binomial_cdf(x, size, p) = pbinon(x, size, p) = P[X <= x]
>
> and P[15 <= X <= 25] = P[X <= 25] - P[X < 15] = P[X <= 25] - P[X <= 14].

You can always try this stuff in Stan, too, using the generated
quantities block (or transformed parameters, though generated
quantities are more efficient).

> In R:
>
> > p <- 0
> > for (i in 15:25) p <- p + choose(50,i) * 0.4^i * (1-0.4)^(50-i)
> > p
> [1] 0.8887012
> > pbinom(25, 50, 0.4) - pbinom(14, 50, 0.4)
> [1] 0.8887012

- Bob

[Martin Šmíra]

Ok, my bad. Sergio was right all the time. We made a mistake in Jags the Jags model to which I was comparing the results.  

Thanks for your patience and explanation.

Martin