RStan significant digits in print()

[Bob Carpenter]

Here's what the default print() looks like for
me in RStan 1.0.2 (latest R, latest Mac OS X):

> print(fit)
Inference for Stan model: dawid-skene.
4 chains: each with iter=500; warmup=250; thin=1; 500 iterations saved.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
pi[1] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 747 1
pi[2] 0.6 0.0 0.0 0.6 0.6 0.6 0.6 0.6 646 1
...

I think the default should have one more
significant digit. The Rhat all print as "1", but
they show up colored on the plot(), so I'm assuming
some of them are actually above 1.1, but you can't
see this. Also, we have enough data for pi here
(25K points for a 7-way multinomial parameter) that I'm
thinking we could use more information on both the values
and the sd and se, which are not actually 0.

I know I can use (digits=xxx), but the woman I was
working with at CCLS didn't know that and just found
the print perplexing as given.

- Bob

[Jiqiang Guo]

The Rhats in the plot are all colored but the intensities are different.  I guess the rhats in the printout is smaller than 1.1.  Otherwise, it would at least be printed as 1.1. 

I do not know what is the issue with the pi here, but their 2.5% and 97.5 percentiles are the same, which indicates the sd is close to 0.  And the more sample sizes we have, the se_mean should be closer to 0. 

--

Jiqiang 

- Bob

--

[Ben Goodrich]

We might want to call format.data.frame() inside the print method for a stanfit object. That does the formatting column-by-column, so if it has to use scientific notation or something for one column, it can use regular formatting for the others.

Ben

[Bob Carpenter]

Default (equiv to digits=1)

> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
> pi[1] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 747 1
> pi[2] 0.6 0.0 0.0 0.6 0.6 0.6 0.6 0.6 646 1

For comparison, here's digits=2, which is actually what I want to see in this example; you
can see that pi[1] and pi[2] are actually getting sampled and that not all Rhat=1.000000...:

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

pi[1] 0.01 0.00 0.00 0.00 0.01 0.01 0.01 0.02 812 1.00
pi[2] 0.61 0.00 0.02 0.58 0.60 0.61 0.62 0.64 487 1.00
pi[3] 0.14 0.00 0.01 0.12 0.13 0.14 0.15 0.17 703 1.00
pi[4] 0.14 0.00 0.01 0.12 0.13 0.14 0.14 0.16 622 1.00
pi[5] 0.04 0.00 0.01 0.03 0.04 0.04 0.05 0.06 417 1.01
pi[6] 0.06 0.00 0.01 0.05 0.06 0.06 0.07 0.08 569 1.00
pi[7] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 892 1.00
theta[1,1,1] 0.69 0.00 0.07 0.56 0.64 0.69 0.74 0.82 795 1.00

Here's digits = 5. It's a bit puzzling that we have Rhat < 1; hopefully it's just floating point shenanigans.

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

pi[1] 0.00902 0.00011 0.00309 0.00427 0.00666 0.00862 0.01099 0.01594 812 0.99985
pi[2] 0.60723 0.00068 0.01508 0.57633 0.59820 0.60667 0.61734 0.63817 487 1.00089
pi[3] 0.14112 0.00045 0.01202 0.12062 0.13211 0.14032 0.14919 0.16548 703 1.00006
pi[4] 0.13718 0.00047 0.01161 0.11592 0.12924 0.13655 0.14445 0.16004 622 1.00046
pi[5] 0.04141 0.00034 0.00693 0.02832 0.03691 0.04085 0.04536 0.05669 417 1.01089
pi[6] 0.06203 0.00030 0.00718 0.04846 0.05717 0.06153 0.06679 0.07716 569 1.00445
pi[7] 0.00202 0.00005 0.00143 0.00020 0.00095 0.00171 0.00280 0.00531 892 0.99901
theta[1,1,1] 0.69148 0.00246 0.06927 0.55581 0.64452 0.69231 0.74228 0.82143 795 0.99852
theta[1,1,2] 0.04851 0.00104 0.03105 0.00616 0.02444 0.04435 0.06540 0.12391 883 1.00023
theta[1,1,3] 0.05009 0.00119 0.03447 0.00681 0.02431 0.04250 0.06898 0.13813 845 1.00007
theta[1,1,4] 0.05241 0.00133 0.03593 0.00439 0.02539 0.04296 0.07444 0.13965 735 0.99977
theta[1,1,5] 0.05287 0.00126 0.03512 0.00572 0.02709 0.04587 0.07053 0.14347 782 0.99672

> The Rhats in the plot are all colored but the intensities are different. I guess the rhats in the printout is
> smaller than 1.1. Otherwise, it would at least be printed as 1.1.

Does R just round down? I'd have thought 1.051 would be printed as 1.1, too.

- Bob

[Jiqiang Guo]

No, in R, I think round(1.051) = 1.1.  There are some  

details in the help of round in R such as

     Note that for rounding off a 5, the IEC 60559 standard is expected

     to be used, ‘_go to the even digit_’.  Therefore ‘round(0.5)’ is

     ‘0’ and ‘round(-1.5)’ is ‘-2’.  However, this is dependent on OS

     services and on representation error (since e.g. ‘0.15’ is not

     represented exactly, the rounding rule applies to the represented

     number and not to the printed number, and so ‘round(0.15, 1)’

     could be either ‘0.1’ or ‘0.2’).

But I do not find that is contradicting to what I said though I do not understand 'go to the even digit'. 

Rhat could be a little less than 1 sometimes; it can happen, not floating point shenanigans, I think. 

--

Jiqiang 

[Bob Carpenter]

On 10/26/12 3:02 PM, Jiqiang Guo wrote:

> Rhat could be a little less than 1 sometimes; it can happen, not floating point shenanigans, I think.

Is this because of the var+ estimator?

I'm pretty sure the "traditional" Gelman and Rubin
R-hat was bounded below at 1.

> No, in R, I think round(1.051) = 1.1.

That's what I'd have expected. I just thought you might've

meant floor() rather than round() when you said:

> I guess the rhats in the printout is
> smaller than 1.1. Otherwise, it would at least be printed as 1.1.

So if I understand correctly, if the rhat are all smaller than 1.05,
they'll all print as 1. And if one has value 1.06, the column will print
values as 1.0 or 1.1? (I don't really care what happens at exactly 1.05.)

- Bob

[Jiqiang Guo]

On Fri, Oct 26, 2012 at 3:26 PM, Bob Carpenter <ca...@alias-i.com> wrote:

On 10/26/12 3:02 PM, Jiqiang Guo wrote:

> Rhat could be a little less than 1 sometimes; it can happen, not floating point shenanigans, I think.

Is this because of the var+ estimator?

I'm pretty sure the "traditional" Gelman and Rubin
R-hat was bounded below at 1.

Someone correct me if I was wrong.  I do not think rhat is always no less than 1. Mathematically, the traditional one is equivalent to what we use now.  I did not find any math inequality that would make sure that B >= M, where B and W are the between- and with-sequence variance.  So I just illustrate by some examples using coda in R. 

> library(coda)

> gelman.diag(mcmc.list(lapply(1:3, function(x) as.mcmc(rnorm(100)))))

Potential scale reduction factors:

     Point est. Upper C.I.

[1,]       1.01       1.02

> gelman.diag(mcmc.list(lapply(1:3, function(x) as.mcmc(rnorm(100)))))

Potential scale reduction factors:

     Point est. Upper C.I.

[1,]          1       1.01

> gelman.diag(mcmc.list(lapply(1:3, function(x) as.mcmc(rnorm(100)))))

Potential scale reduction factors:

     Point est. Upper C.I.

[1,]       1.06       1.21

> gelman.diag(mcmc.list(lapply(1:3, function(x) as.mcmc(rnorm(100)))))

Potential scale reduction factors:

     Point est. Upper C.I.

[1,]      0.992      0.995

> gelman.diag(mcmc.list(lapply(1:3, function(x) as.mcmc(rnorm(100)))))

Potential scale reduction factors:

     Point est. Upper C.I.

[1,]          1       1.02

> gelman.diag(mcmc.list(lapply(1:3, function(x) as.mcmc(rnorm(100)))))

Potential scale reduction factors:

     Point est. Upper C.I.

[1,]      0.998       1.02

 

No, in R, I think round(1.051) = 1.1.

That's what I'd have expected.  I just thought you might've

meant floor() rather than round() when you said:

> I guess the rhats in the printout is
> smaller than 1.1.  Otherwise, it would at least be printed as 1.1.

So if I understand correctly, if the rhat are all smaller than 1.05,
they'll all print as 1.  And if one has value 1.06, the column will print
values as 1.0 or 1.1?  (I don't really care what happens at exactly 1.05.)

It should be 1.1 as the following shows

> print(round(1.06, 1))

[1] 1.1

 

- Bob

--

[Bob Carpenter]

On 10/26/12 7:19 PM, Jiqiang Guo wrote:

> Someone correct me if I was wrong. I do not think rhat is always no less than 1. Mathematically, the traditional one is
> equivalent to what we use now. I did not find any math inequality that would make sure that B >= M, where B and W are
> the between- and with-sequence variance. So I just illustrate by some examples using coda in R.

You're right (I'm now feeling like a guy on Jeopardy
who keeps buzzing in with the wrong answer)!

I'd never seen values of r-hat below 1 before, but you
wouldn't expect to with diffuse starting points in MCMC.

So I came up with a simple example:

withins:

> sd(c(1,1,2))
[1] 0.5773503

> sd(c(1,2,2))
[1] 0.5773503

between:

> sd(c(1,1,1,2,2,2))
[1] 0.5477226

The n/(n-1) in the sd calc matters (the within and between
are different sizes). With more samples, it flips back to
what you'd expect:

> sd(c(1,1,1,1,2,2))
[1] 0.5163978

> sd(c(1,1,1,1,1,1,2,2,2,2,2,2))
[1] 0.522233


Andrew and Kenny say this in their MCMC Handbook paper:

> At convergence, the chains will have mixed, so that the distribution of the simulations between and within chains
> will be identical, and the ratio R.hat should equal 1. If R.hat is greater than 1, this implies that the chains have
> not fully mixed and that further sim- ulation might increase the precision of inferences.

With r-hat < 1, the within-chain estimates don't look
like the between-chain estimate.

> mean(c(1,1,2))
[1] 1.333333

> mean(c(2,2,1))
[1] 1.666667

> mean(c(1,1,1,2,2,2))
[1] 1.5

So do we want to say when |rhat - 1| > epsilon we haven't converged?

Not that I expect this to come up much.

- Bob