"nan" values for Rhat in pystan

[Jesse Bloom]

I sometimes get Rhat values of "nan" when using PyStan. This happens for certain data sets repeatedly. I think I get these values for variables that are actually very close to zero and seem well converged based on the effective sample size (n_eff) values. Is this a bug in stan?

[Bob Carpenter]

If R-hat is NaN and n_eff is not NaN, it's almost certainly
a precision problemin the arithmetic.

Could you send us a model and data set for which it happens so
we can diagnose?

There may be loss of arithmetic precision in R-hat compared to
what we could do ideally. The problem is that floating point
subtraction of values close to each other results in so-called
catastrophic cancellation, losing an arbitrary amount of precision.
Something like Welford's algorithm is required to maintain stability
for sample variance calculations, and I'm not sure how PyStan has implemented
R-hat.

For example, this model is problematic arithmetically:

fit <- stan(model_code = "parameters { real y; } model { y ~ normal(1,1e-20); }")

and produces

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
y 1 0.0 0 1 1 1 1 1 4000 Inf
lp__ -15407440 422002.9 26689805 -61629758 -15407440 0 0 0 4000 1.500106e+14


There's really nothing we can do about such models numerically. What
you have to do from a user perspective is use a parameterizations

z = y - 1,

which can be handled stably because you can represent numbers close to zero well.

parameters { real z; } model { z ~ normal(0,1e-20); }

The problem is that there's still no way to represent y = z + 1, but
you can get close to log(y) by using

log_y <- log1p(z);

We're gradually pulling all these implementations back into
the core Stan C++ code, so at least we'll only have one fix to make
for issues like this in the near future.

- Bob

P.S. I'm going to drop the above discussion into a new section on
arithmetic in the problematic posteriors chapter of the manual
or maybe an entire chapter on arithmetic stability of various operators.

> On Dec 18, 2014, at 11:40 AM, Jesse Bloom <jesse...@gmail.com> wrote:
>
> I sometimes get Rhat values of "nan" when using PyStan. This happens for certain data sets repeatedly. I think I get these values for variables that are actually very close to zero and seem well converged based on the effective sample size (n_eff) values. Is this a bug in stan?
>

> --
> 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/d/optout.

[Allen B. Riddell]

Hi Jesse,

Yes, this is an issue. While R will let you divide by zero (or nearly
zero) quickly and without any problem (other than throwing up an 'Inf'),
Python raises an exception. The current way PyStan handles this probably
isn't the best. If you have a simple case where nan isn't appropriate,
would you mind opening an issue on PyStan's github site?

On the other hand, if there really is zero variance among the chains
with respect to some parameter, rhat and n_eff should be 'nan'.

Thanks,

-a

[Bob Carpenter]

Are you using double-precision floating point? As far as I
know, R's just using standard IEEE floating point.

The usual problems are with subtraction, not division, unless
you get so close to zero that the division overflows to +inf.
But if that happens for two numbers, it should happen in R, too.

And then once you get inf, everything is inf or NaN. You can't
divide inf/inf or subract (inf - inf), although you can add (inf + inf)
and get inf as long as the infs are the right sign.

- Bob

[Allen B. Riddell]

Hi Bob,

Yeah, Python throws an error we have to catch:

> In [4]: 1 / 1e-500
> ---------------------------------------------------------------------------
> ZeroDivisionError Traceback (most recent call last)
> <ipython-input-4-ba7ce1d0dee4> in <module>()
> ----> 1 1 / 1e-500
>
> ZeroDivisionError: float division by zero

rhat and n_eff are bottlenecks on huge models so there's been a lot of
churn trying to get it right AND fast. Originally I had everything going
through Stan and C++ but handing things off to C++ and then handing it
back to Python was too inefficient. If Stan could expose some way for
RStan and PyStan to calculate rhat and n_eff, that would be fantastic.
Right now I believe there are 3 implementations of the rhat and n_eff
calculations (one for each language).

I'd definitely welcome help on this, if anyone is interested.

Thanks,

-a

[Paula Sarkis]

Hello, 

I am having the same problem. The trace plot looks perfect however the value of Rhat is nan. 

Can you please give me some details on how to evaluate n_eff for each parameter. Also, how to evaluate lp. I've been having the 'nan' problem with rhat many times now, I want to double check whether my model is ill defined or if the problem is from python. 

Cheers,

Paula

[Bob Carpenter]

The manual provides the formula in the MCMC chapter.

I don't know that there's a way for R-hat to be NaN without
n_eff being NaN.

If PyStan lets you dump the output in the CSV format which
CmdStan can read, you could check that way.

- Bob

> To post to this group, send email to stan-...@googlegroups.com.

[Paula Sarkis]

I have around 100 parameters. All the parameters have n_eff = 6668 and Rhat = 1.0 except for only 1 parameter (the most important one). This parameter has n_eff = 6668 but Rhat = nan

I am running MCMC with 10,000 iterations, warmup = 5,000 and with 4 chains. So that's a total of 20,000 draws. n_eff of lp is 3938 and Rhat = 1.0.

Any idea why I am getting Rhat  = nan?

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/hn4W_p8j3fs/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Allen B. Riddell]

Dear Paula,

PyStan's rhat calculation has some bugs associated with it. In the
future, PyStan will use the same rhat code as CmdStan uses.

There are known issues when the variance of the parameter is zero but it
sounds like that's not what's happening here. If you can reproduce the
problem with model code and data you're willing to share, I can look
into it. If that's not possible, I'd suggest trying to run the model in
RStan and seeing if rhat is displayed.

Best,

Allen

[Paula Sarkis]

Dear Allen, 

I appreciate your help, thank you!

Please find attached the data and the model. 

Here's a quick description of the model:

I want to fit the data points with a line. The data has measurement errors. So xerr and yerr is simply the uncertainties on the parameters.

Cheers!

Paula

[Paula Sarkis]

Dear Allen, 

I would like to note that when I run the model on a 'toy model' the results match the true values of the parameters. Also, I get rhat = 1.0 for all the parameters. It is only when I run the model on my data I get rhat = nan for 1 parameter (which is the most important one) and for the other parameters rhat = 1.0.

Best, 

Paula

[Jason Tilley]

I am also having this issue for some of my parameters in PyStan. Rhat = nan, but n_eff=iterations. My model code is:

import numpy as np

import matplotlib.pyplot as plt

import pystan

# model predicting weight and height from waist size and shoe size.
# fake data

N = 100 # number per group

# men

height_m = np.random.normal(5.7, 0.3, N)

shoe_size_m = height_m * 3.5 + np.random.normal(0,0.5,N) - 10

weight_m = height_m ** 2.8 + (np.random.gamma(3,2,N) - 2 ** 2) * 18

waist_size_m = weight_m ** 0.33 * height_m + 4 + np.random.normal(0,0.5,N)

is_male_m = np.ones(N)

# women

height_w = np.random.normal(5.2, 0.3, N)

shoe_size_w = height_w * 3.4 + np.random.normal(0,0.5,N) - 10

weight_w = height_w ** 2.8 + (np.random.gamma(3,2,N) - 2 ** 2) * 10 + 20

waist_size_w = weight_w ** 0.33 * height_w + np.random.normal(0,0.7,N)

is_male_w = np.ones(N) + 1

height = np.append(height_m, height_w)

shoe_size = np.append(shoe_size_m, shoe_size_w)

weight = np.append(weight_m, weight_w)

waist_size = np.append(waist_size_m, waist_size_w)

group = np.append(is_male_m, is_male_w).astype(int)

y = np.matrix([list(height), list(np.log(weight))]).T

x = np.matrix([list(waist_size), list(shoe_size)]).T

model_dat = ({'N':len(y),

             'J': np.shape(x)[1],

             'K': np.shape(y)[1],

             'L': len(np.unique(group)),

             'y': y,

             'x': x,

             'random': group,

             })

model_code = """

    data {

      int<lower=0> N; //rows of data

      int<lower=1> J; // number of covariates

      int<lower=1> K; // number of responses

      int<lower=1> L; // number of groups

      vector[J] x[N]; // design matrix (predictors)

      vector[K] y[N]; // response (outcomes)

      int<lower=1, upper=L> random[N];

    }

    parameters {

      vector[K] intercept; // Population intercept

      matrix[K,J] beta; // Population slope parameters

      cholesky_factor_corr[K] omega;

      cholesky_factor_corr[L] omega_u;

      vector<lower=0>[K] sigma; // population standard deviation

      vector[L] u; // random intercepts

      vector<lower=0>[L] sigma_u; // random intercepts SD

    }

    model {

      vector[K] mu[N];

      matrix[K,K] Sigma;

      matrix[L,L] Sigma_u;

      vector[L] zeros;

      

      Sigma <- diag_pre_multiply(sigma, omega);

      Sigma_u <- diag_pre_multiply(sigma_u, omega_u);

      

      for (l in 1:L)

        zeros[l] <- 0;

      

      u ~ multi_normal_cholesky(zeros, Sigma_u);

        

      for (n in 1:N)

        mu[n] <- beta * x[n] + intercept + u[random[N]];

        

      // Priors

      to_vector(beta) ~ normal(0, 5);

      omega ~ lkj_corr_cholesky(4);

      omega_u ~ lkj_corr_cholesky(4);

      sigma ~ cauchy(0, 2.5);

      sigma_u ~ cauchy(0, 2.5);

      

      // Data Model

      y ~ multi_normal_cholesky(mu, Sigma);

    }

    """

fit = pystan.stan(model_code=model_code, data=model_dat,

                  iter=1000, chains=2,

                  model_name='height_weight',

                  sample_file='/Users/username/Desktop/sample')

My results are:

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

intercept[0]   6.57    2.76   10.7   -8.2   2.41   3.61   6.47  41.43   15.0   1.15

intercept[1]   7.01    2.76   10.7  -7.76   2.82   4.06   6.94  41.86   15.0   1.15

beta[0,0]      0.01  2.0e-4 2.9e-3 7.4e-3   0.01   0.01   0.02   0.02  223.0   1.01

beta[1,0]      0.08  2.2e-4 3.4e-3   0.07   0.07   0.08   0.08   0.08  231.0    1.0

beta[0,1]      0.19  6.3e-4 9.4e-3   0.17   0.18   0.19    0.2   0.21  223.0   1.01

beta[1,1]     -0.14  7.7e-4   0.01  -0.16  -0.14  -0.14  -0.13  -0.11  223.0    1.0

omega[0,0]      1.0     0.0    0.0    1.0    1.0    1.0    1.0    1.0 1000.0    nan

omega[1,0]     -0.2  4.5e-3   0.07  -0.33  -0.25   -0.2  -0.15  -0.07  231.0    1.0

omega[0,1]      0.0     0.0    0.0    0.0    0.0    0.0    0.0    0.0 1000.0    nan

omega[1,1]     0.98  9.4e-4   0.01   0.94   0.97   0.98   0.99    1.0  227.0    1.0

omega_u[0,0]    1.0     0.0    0.0    1.0    1.0    1.0    1.0    1.0 1000.0    nan

omega_u[1,0]  -0.01    0.02   0.33  -0.64  -0.26  -0.01   0.24    0.6  252.0    1.0

omega_u[0,1]    0.0     0.0    0.0    0.0    0.0    0.0    0.0    0.0 1000.0    nan

omega_u[1,1]   0.94  4.4e-3   0.07   0.73   0.92   0.97   0.99    1.0  278.0    1.0

sigma[0]       0.13  5.0e-4 7.1e-3   0.12   0.13   0.13   0.14   0.15  200.0    1.0

sigma[1]       0.16  5.2e-4 8.3e-3   0.15   0.16   0.16   0.17   0.18  252.0   1.01

u[0]          -0.59    1.23   8.52 -25.04  -1.77   0.01   1.65  18.21   48.0    1.0

u[1]          -3.19    2.76   10.7 -38.04  -3.08  -0.22    1.0  11.58   15.0   1.15

sigma_u[0]     7.01    1.53  13.79   0.16   1.26   2.88   6.59  43.86   81.0   1.04

sigma_u[1]     7.46    2.01   12.4    0.1   1.11   2.71   7.84  44.05   38.0   1.04

lp__         565.99    0.47   3.33 558.57 563.81  566.2 568.46 571.71   51.0   1.05

Not great. Also, the fit.plot() function throws this error:

>>> fit.plot()

Traceback (most recent call last):

  File "<stdin>", line 1, in <module>

  File "stanfit4height_weight_2bb6f32b402d0b58a5f60e9816a2700b_3413209832667885730.pyx", line 491, in stanfit4height_weight_2bb6f32b402d0b58a5f60e9816a2700b_3413209832667885730.StanFit4Model.plot (/var/folders/b1/zqkx5fbx7z933zhdmdbxn3s00000gn/T/tmpQ06aBv/stanfit4height_weight_2bb6f32b402d0b58a5f60e9816a2700b_3413209832667885730.cpp:9487)

  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/pystan/plots.py", line 19, in traceplot

    return plots.traceplot(fit.extract(), pars, **kwargs)

  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/pystan/external/pymc/plots.py", line 41, in traceplot

    kdeplot_op(ax[i, 0], d)

  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/pystan/external/pymc/plots.py", line 60, in kdeplot_op

    ax.plot(x, density(x))

  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/stats/kde.py", line 204, in evaluate

    raise ValueError(msg)

ValueError: points have dimension 1, dataset has dimension 1000

Any ideas?

[Ben Goodrich]

On Friday, January 22, 2016 at 7:35:47 PM UTC-5, Jason Tilley wrote:

I am also having this issue for some of my parameters in PyStan. Rhat = nan, but n_eff=iterations.

For a Cholesky factor, L, the L[i,j] elements for j > i are fixed at 0. Moreover, for a Cholesky factor of a correlation matrix, the top left element is fixed at 1. Therefore, no function of their "distributions" --- including Rhat --- is well-defined.

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

omega[0,0]      1.0     0.0    0.0    1.0    1.0    1.0    1.0    1.0 1000.0    nan

omega[1,0]     -0.2  4.5e-3   0.07  -0.33  -0.25   -0.2  -0.15  -0.07  231.0    1.0

omega[0,1]      0.0     0.0    0.0    0.0    0.0    0.0    0.0    0.0 1000.0    nan

omega[1,1]     0.98  9.4e-4   0.01   0.94   0.97   0.98   0.99    1.0  227.0    1.0

omega_u[0,0]    1.0     0.0    0.0    1.0    1.0    1.0    1.0    1.0 1000.0    nan

omega_u[1,0]  -0.01    0.02   0.33  -0.64  -0.26  -0.01   0.24    0.6  252.0    1.0

omega_u[0,1]    0.0     0.0    0.0    0.0    0.0    0.0    0.0    0.0 1000.0    nan

omega_u[1,1]   0.94  4.4e-3   0.07   0.73   0.92   0.97   0.99    1.0  278.0    1.0

 
Ben

[Ben Goodrich]

On Friday, January 22, 2016 at 7:35:47 PM UTC-5, Jason Tilley wrote:

I am also having this issue for some of my parameters in PyStan. Rhat = nan, but n_eff=iterations.

For a Cholesky factor, L, the L[i,j] elements for j > i are fixed at 0. Moreover, for a Cholesky factor of a correlation matrix, the top left element is fixed at 1. Therefore, no function of their "distributions" --- including Rhat --- is well-defined.

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

omega[0,0]      1.0     0.0    0.0    1.0    1.0    1.0    1.0    1.0 1000.0    nan

omega[1,0]     -0.2  4.5e-3   0.07  -0.33  -0.25   -0.2  -0.15  -0.07  231.0    1.0

omega[0,1]      0.0     0.0    0.0    0.0    0.0    0.0    0.0    0.0 1000.0    nan

omega[1,1]     0.98  9.4e-4   0.01   0.94   0.97   0.98   0.99    1.0  227.0    1.0

omega_u[0,0]    1.0     0.0    0.0    1.0    1.0    1.0    1.0    1.0 1000.0    nan

omega_u[1,0]  -0.01    0.02   0.33  -0.64  -0.26  -0.01   0.24    0.6  252.0    1.0

omega_u[0,1]    0.0     0.0    0.0    0.0    0.0    0.0    0.0    0.0 1000.0    nan

omega_u[1,1]   0.94  4.4e-3   0.07   0.73   0.92   0.97   0.99    1.0  278.0    1.0

 
Ben

[Jason Tilley]

Thanks. So I assume this is ok then? I get better Rhat with a bigger sample size, so that is good.

[Ben Goodrich]

On Friday, January 22, 2016 at 9:35:45 PM UTC-5, Jason Tilley wrote:

Thanks. So I assume this is ok then? I get better Rhat with a bigger sample size, so that is good.

For those "margins" of the posterior distribution, yes, but only because they aren't margins of the posterior distribution.

[Bob Carpenter]

They're not margins in the sense that they're deterministic.
I just figured it'd be easier to scan that way, but it's
definitely wasteful.

Right now, the problem in printing is that there are two
ways things can have the same value in every draw. One is
to be a discrete generated quantity that never changes. This
can happen in event probability calculations if the probability
is very close to 0 or 1 or the chains aren't run for many
iterations (you need many more iterations to estimate tail
probabilities well).

So I don't know if there's a good way to tell the ones that
are constant by design (diagonal of a correlation matrix) and
ones that are random quantities, but constant in a chain
because every draw has produced the same value.

- Bob