Neff, Rhat for constant looks strange

tau_k[1]         1.000000e+00 0.0000000000  0.00000000  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00 1000.000000       NaN
tau_k[2]         1.137104e+00 0.0038440625  0.06419562  1.015457e+00  1.093082e+00  1.137091e+00  1.181873e+00  1.258226e+00  278.888019 1.0194419
tau_k[3]         9.750000e-01 0.0000000000  0.00000000  9.750000e-01  9.750000e-01  9.750000e-01  9.750000e-01  9.750000e-01    1.003009 0.9979980

vector make_tau(real r) {
vector[3] x;
x[1] = 1.0;
x[2] = r;
x[3] = 0.975;
return x;
}

parameter

real<lower=0> tau_k_raw;

transformed parameter:

vector<lower=0>[3] tau_k = make_tau(tau_k_raw);

tau_k_raw ~ normal(0, 1);

neff, Rhat for tau_k[3] looks strange. No question, I can get around it. I don’t need a
solution. Although it may be an indicator that something internal needs some checking.

options(max.print=1000000, width = 999)
sink(paste(fname, “.diagnosis”, sep=“”))
looic ← loo(extract_log_lik(lm_waic))
print(looic)
print(summary(lm_waic))
sink()

Loading required package: ggplot2
Stackoverflow is a great place to get help:
Newest 'ggplot2' Questions - Stack Overflow.
Loading required package: StanHeaders
rstan (Version 2.16.2, packaged: 2017-07-03 09:24:58 UTC, GitRev: 2e1f913d3ca3)
For execution on a local, multicore CPU with excess RAM we recommend calling
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

Ubuntu 16.04. gcc 5.4.0. default installation, no tweaks.

Andre

1 Like

Yeah that’s true. yhat[1] and yhat[3] are both constants, but one produces a NaN Rhat and one produces an Rhat < 1.0. And the n_effs are different as well, even though they are both constants close to each other.

The output of the first seems reasonable, but the second is weird. I do not know what’s going on here.

0.975 can’t be represented exactly so there is some numerical noise to how it gets stored and transmitted back to the interface.

But wouldn’t it be the same rounding every time?

Apparently not.

Touche.

Maybe fft - ifft plays a trick for neff.

And var for rhat? Are there lower limits for variances in RStan / CmdStan like PyStan has?

The estimated Rhat can be < 1, even with random variables.

Line 211

if var_within < pystan.constants.EPSILON:
    srhat = float('nan')

I never managed to get the calculation for rhat and neff in PyStan to match the one in RStan exactly. It would be very, very nice to have these calculations done by Stan in C++ rather than in the interfaces. It bothers me that the results aren’t the same across interfaces.

IIRC some of this code addresses a difference between R and Python: R doesn’t complain about dividing a float by zero whereas Python raises a ZeroDivision exception.

This seems to be the closest related post I could find (also possibly related: Output for quantities of interest that are actually constants). I have a general question about posterior samples of constants. Sometimes they seem to yield NaN neffs and Rhats, and in the same model other constants has small (but defined) neff and Rhat of 1. The model converges fine and so this is not a critical issue (I assume). I’m just curious what brings about the consistent output.

E.g., in the following output for which m_0 and S_ are user-provided means and covariance matrices of bivariate Gaussians (for two categories each, indicated by the first array dimension), and bias is a user-provided simplex of prior probabilities (for the two categories):

...
m_0[1,1]                 14.00     NaN  0.00    14.00    14.00    14.00    14.00    14.00   NaN  NaN
m_0[1,2]                138.35    0.00  0.00   138.35   138.35   138.35   138.35   138.35     2    1
m_0[2,1]                 61.00     NaN  0.00    61.00    61.00    61.00    61.00    61.00   NaN  NaN
m_0[2,2]                157.91    0.00  0.00   157.91   157.91   157.91   157.91   157.91     2    1
S_0[1,1,1]               81.00     NaN  0.00    81.00    81.00    81.00    81.00    81.00   NaN  NaN
S_0[1,1,2]              106.89    0.00  0.00   106.89   106.89   106.89   106.89   106.89     2    1
S_0[1,2,1]              106.89    0.00  0.00   106.89   106.89   106.89   106.89   106.89     2    1
S_0[1,2,2]            11658.50    0.00  0.00 11658.50 11658.50 11658.50 11658.50 11658.50     2    1
S_0[2,1,1]             6561.00     NaN  0.00  6561.00  6561.00  6561.00  6561.00  6561.00   NaN  NaN
S_0[2,1,2]             -291.27    0.00  0.00  -291.27  -291.27  -291.27  -291.27  -291.27     2    1
S_0[2,2,1]             -291.27    0.00  0.00  -291.27  -291.27  -291.27  -291.27  -291.27     2    1
S_0[2,2,2]              399.11    0.00  0.00   399.11   399.11   399.11   399.11   399.11     2    1
bias[1]                   0.50     NaN  0.00     0.50     0.50     0.50     0.50     0.50   NaN  NaN
bias[2]                   0.50     NaN  0.00     0.50     0.50     0.50     0.50     0.50   NaN  NaN
...

I’m using rstan 2.21.2 (but this output inconsistency has been the case as long as I can remember). I wonder whether it is something that rstan might be able to catch? Thank you and sorry for reviving this old topic.

platform       x86_64-apple-darwin17.0     
arch           x86_64                      
os             darwin17.0                  
system         x86_64, darwin17.0          
status                                     
major          4                           
minor          0.2                         
year           2020                        
month          06                          
day            22                          
svn rev        78730                       
language       R                           
version.string R version 4.0.2 (2020-06-22)
nickname       Taking Off Again

Can you provide Rdata file of the values which produce 1? Test also you get the same result after loading the values from the file you are going to send. I suspect floating point accuracy problem as was found that Eigen’s acov estimate can produce.

I assume you mean a file with the posterior samples? It seems I can’t upload RData files here, but here’s a link to a dropbox file:

https://www.dropbox.com/s/5vb2wjorq5zoo41/samples.RData?dl=0 (includes object x which is a tidybayes::spread_draws() of the m_0 and S_0 samples along the variables cat, cue1, cue2)

and in case you mean the model itself:

There is zero variance among those values:

> library(tidyverse)
>  x %>% filter(cat == 1, cue1 == 1) %>% summarise(m_0 = var(m_0))
    `summarise()` regrouping output by 'cat', 'cue1' (override with `.groups` argument)
    # A tibble: 2 x 4
    # Groups:   cat, cue1 [1]
        cat  cue1 cue2    m_0
      <int> <int> <fct> <dbl>
    1     1     1 VOT       0
    2     1     1 f0        0

The R-hats are all slightly below 1 (0.9996666) before and after loading. Apologies if this doesn’t address your suggestion (and thank you for your reply).

1 Like

The minimal set so that I can reproduce what you observe.

I haven’t used tidybayes. How do you compute Rhat and ESS for this object? is the output in your earlier post from tidybayes? Do you know which code it does use for Rhat and ESS computation? Its own or code from RStan or code from posterior package?

If you load the model object (from the second RData file I linked), and you summarize it, it’ll give you rhat, Neff, etc. You can also use rhat(model) once you have loaded it.

I only brought up tidybayes as one way to extract the posterior samples, and to confirm that they have zero variance (as they should, since these are variables that are constants).

I get

> summary(ibbu.exp2.with_prior_no_lapse)
          Length            Class             Mode 
               1 NIW_ibbu_stanfit               S4 

It seems I need a bit more info how to replicate your result,

Looks like the class got mangled. You can fix it with:

library(rstan)

load("model.RData")

new_stanfit_obj <- ibbu.exp2.with_prior_no_lapse
class(new_stanfit_obj) <- "stanfit"
new_stanfit_obj

I think what @tiflo is referring to can be seen with:

library(abind)
acorn(extract(new_stanfit_obj, "m_0", permuted = FALSE), n = 5, r = 3)

which makes

seem likely.

1 Like

Correct, @hhau (and thank you). This is a class that adds a few slots to the stanfit class, but these slots are independent of the issue I raised.

@avehtari, I didn’t think of this issue (since I have the library that creates the new class, it didn’t show on my side). But the issue with the Neff and rhats should still show the way that hhau suggested. Does it for you? If not I can recreate and example with just rstan.

Thanks @hhau. I can now test different codes. Pinging @bgoodri

With RStan 2.21.2

print shows the inconsistent behavior

print(new_stanfit_obj)
...
m_0[1,1]                 14.00     NaN  0.00    14.00    14.00    14.00    14.00    14.00   NaN  NaN
m_0[1,2]                138.35    0.00  0.00   138.35   138.35   138.35   138.35   138.35     2    1
m_0[2,1]                 61.00     NaN  0.00    61.00    61.00    61.00    61.00    61.00   NaN  NaN
m_0[2,2]                157.91    0.00  0.00   157.91   157.91   157.91   157.91   157.91     2    1
...

rstan::monitor shows consistently Rhat=1 and ESSs are equal to the sample size (interestingly the row order is different).

rstan::monitor(new_stanfit_obj)
...
m_0[1,1]                 14.0    14.0    14.0    14.0  0.0     1    12000    12000
m_0[2,1]                 61.0    61.0    61.0    61.0  0.0     1    12000    12000
m_0[1,2]                138.3   138.3   138.3   138.3  0.0     1    12000    12000
m_0[2,2]                157.9   157.9   157.9   157.9  0.0     1    12000    12000
...

Rhat=1 and ESS equal to the sample size would be fine if we know these are constants. However, what Stan returns doesn’t contain information which variables are truly constant and which are not truly constant. It’s possible that for truly non-constant we observe just by chance that all posterior draws are equal. In that case Rhat and ESS are not well defined. As we can’t make difference between these cases it was decided that for all equal draws Rhat and ESS should be NA. New posterior package does this

posterior::summarise_draws(posterior::as_draws(new_stanfit_obj))
...
 6 m_0[1,1]          14       14       0         0          14       14     NA         NA       NA 
 7 m_0[2,1]          61       61       0         0          61       61     NA         NA       NA 
 8 m_0[1,2]         138.     138.      0         0         138.     138.    NA         NA       NA 
 9 m_0[2,2]         158.     158.      0         0         158.     158.    NA         NA       NA 
...

CmdStanR is already using the posterior package. RStan probably will start using it at some point.

1 Like