# Neff, Rhat for constant looks strange

**URL:** <https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011>\
**Category:** Modeling\
**Created:** [September 26, 2017, 2:17am UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011 "2017-09-26T02:17:23Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Andre\_Pfeuffer](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andre_pfeuffer/32/4939_2.png) [@Andre\_Pfeuffer](https://discourse.mc-stan.org/u/Andre_Pfeuffer)\
**Post date:** [September 26, 2017, 2:17am UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/1 "2017-09-26T02:17:23Z")

</div>

```
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](http://stackoverflow.com/tags/ggplot2).  
> 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

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [September 27, 2017, 12:50am UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/2 "2017-09-27T00:50:40Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [September 27, 2017, 1:08am UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/3 "2017-09-27T01:08:29Z")

</div>

> [@bbbales2](#):
>
> 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.

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [September 27, 2017, 1:09am UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/4 "2017-09-27T01:09:25Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [September 27, 2017, 1:09am UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/5 "2017-09-27T01:09:48Z")

</div>

Apparently not.

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [September 27, 2017, 1:11am UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/6 "2017-09-27T01:11:18Z")

</div>

Touche.

---

<div class="post-metadata">

**Author:** ![ahartikainen](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ahartikainen/32/321_2.png) [@ahartikainen](https://discourse.mc-stan.org/u/ahartikainen)\
**Post date:** [September 27, 2017, 7:01pm UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/7 "2017-09-27T19:01:29Z")

</div>

Maybe fft - ifft plays a trick for neff.

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

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [September 27, 2017, 7:02pm UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/8 "2017-09-27T19:02:57Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ahartikainen](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ahartikainen/32/321_2.png) [@ahartikainen](https://discourse.mc-stan.org/u/ahartikainen)\
**Post date:** [September 27, 2017, 7:09pm UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/9 "2017-09-27T19:09:24Z")

</div>

Line 211

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

```

> <https://github.com/stan-dev/pystan/blob/develop/pystan/_chains.pyx>

---

<div class="post-metadata">

**Author:** ![ariddell](https://avatars.discourse-cdn.com/v4/letter/a/57b2e6/32.png) [@ariddell](https://discourse.mc-stan.org/u/ariddell)\
**Post date:** [September 28, 2017, 2:32pm UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/10 "2017-09-28T14:32:16Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ariddell](https://avatars.discourse-cdn.com/v4/letter/a/57b2e6/32.png) [@ariddell](https://discourse.mc-stan.org/u/ariddell)\
**Post date:** [September 28, 2017, 2:33pm UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/11 "2017-09-28T14:33:57Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![tiflo](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/tiflo/32/3044_2.png) [@tiflo](https://discourse.mc-stan.org/u/tiflo)\
**Post date:** [August 16, 2020, 4:18pm UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/12 "2020-08-16T16:18:17Z")

</div>

This seems to be the closest related post I could find (also possibly related: [Output for quantities of interest that are actually constants](https://discourse.mc-stan.org/t/output-for-quantities-of-interest-that-are-actually-constants/7948)). 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
```

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [August 16, 2020, 4:51pm UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/13 "2020-08-16T16:51:19Z")

</div>

> [@tiflo](#):
>
> other constants has small (but defined) neff and Rhat of 1.

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.

---

<div class="post-metadata">

**Author:** ![tiflo](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/tiflo/32/3044_2.png) [@tiflo](https://discourse.mc-stan.org/u/tiflo)\
**Post date:** [August 17, 2020, 12:06am UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/14 "2020-08-17T00:06:57Z")

</div>

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](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:

> **[model.RData](https://www.dropbox.com/s/13rnxxdnz157mop/model.RData?dl=0)**
>
> Shared with Dropbox

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

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [August 18, 2020, 11:57am UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/15 "2020-08-18T11:57:51Z")

</div>

> [@tiflo](#):
>
> I assume you mean a file with the posterior samples?

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

> [@tiflo](#):
>
> [Dropbox - File Deleted - Simplify your life](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)

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?

---

<div class="post-metadata">

**Author:** ![tiflo](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/tiflo/32/3044_2.png) [@tiflo](https://discourse.mc-stan.org/u/tiflo)\
**Post date:** [August 18, 2020, 11:09pm UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/16 "2020-08-18T23:09:23Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [August 25, 2020, 7:39pm UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/17 "2020-08-25T19:39:39Z")

</div>

> [@tiflo](#):
>
> If you load the model object (from the second RData file I linked), and you summarize it, it’ll give you rhat, Neff, etc.

I get

```no-highlight
> 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,

---

<div class="post-metadata">

**Author:** ![hhau](https://avatars.discourse-cdn.com/v4/letter/h/2bfe46/32.png) [@hhau](https://discourse.mc-stan.org/u/hhau)\
**Post date:** [August 26, 2020, 8:16am UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/18 "2020-08-26T08:16:29Z")

</div>

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

```no-highlight
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:

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

```

which makes

> [@avehtari](#):
>
> floating point accuracy problem as was found that Eigen’s acov estimate can produce.

seem likely.

---

<div class="post-metadata">

**Author:** ![tiflo](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/tiflo/32/3044_2.png) [@tiflo](https://discourse.mc-stan.org/u/tiflo)\
**Post date:** [August 26, 2020, 12:01pm UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/19 "2020-08-26T12:01:16Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [August 27, 2020, 2:34pm UTC](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011/20 "2020-08-27T14:34:31Z")

</div>

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.

[Next page](https://discourse.mc-stan.org/t/neff-rhat-for-constant-looks-strange/2011.md?page=2)