# Slightly different modelling changes the fitting

I just wanted to fit a 2nd order polynomial. In the first exaple I am having no problem:

``````import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np

A0=0.5; A1=1.5; A2=-0.2; A3=-0.008; A4=0.00025
def func(xx):
return A0+A1*xx+A2*xx*xx#+A3*(x**3)+A4*(x**4)

x=10*np.random.rand(100); x=np.sort(x)
fx=func(x);
yerror=np.random.rand(len(x))
y=np.random.normal(fx,scale=yerror)

np.random.seed(101)

model = """
data {
int<lower=0> N;
vector[N] x;
vector[N] y;
vector[N] yerror;
int <lower=0> NP;
}
parameters {
real a[NP];
}
model {
for(i in 1:N){
target+=normal_lpdf(y[i]|a + a * x[i] + a * x[i]*x[i], yerror[i]);

}

}
"""

NP=3;
data = {'N': len(x), 'x': x, 'y': y,'yerror':yerror,'NP':NP}
sm = pystan.StanModel(model_code=model)
fit = sm.sampling(data=data, iter=2000, chains=4, warmup=400, thin=3, seed=101)
print(fit)
``````

The output is normal:

``````INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_9e171b5a8aec523e35681175cfb02420 NOW.
Inference for Stan model: anon_model_9e171b5a8aec523e35681175cfb02420.
4 chains, each with iter=2000; warmup=400; thin=3;
post-warmup draws per chain=534, total post-warmup draws=2136.

mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
a   1.52  8.4e-4   0.02   1.48   1.51   1.52   1.54   1.57    877    1.0
a   -0.2  6.6e-5 2.1e-3  -0.21   -0.2   -0.2   -0.2   -0.2    979    1.0
a   0.44  2.3e-3   0.07   0.31   0.39   0.44   0.49   0.57    844    1.0
lp__ -50.18    0.03   1.21 -53.35  -50.7 -49.84 -49.31 -48.79   1385    1.0

Samples were drawn using NUTS at Wed Nov 20 17:18:26 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
``````

But when I try to do the same with declaring a as a new parameter, I am getting surprising warnings and Rhat is not becoming 1.0. See the following code:

``````import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np

A0=0.5; A1=1.5; A2=-0.2; A3=-0.008; A4=0.00025
def func(xx):
return A0+A1*xx+A2*xx*xx#+A3*(x**3)+A4*(x**4)

x=10*np.random.rand(100); x=np.sort(x)
fx=func(x);
yerror=np.random.rand(len(x))
y=np.random.normal(fx,scale=yerror)

np.random.seed(101)

model = """
data {
int<lower=0> N;
vector[N] x;
vector[N] y;
vector[N] yerror;
int <lower=0> NP;
}
parameters {
real shift;
real a[NP];
}
model {
for(i in 1:N){
target+=normal_lpdf(y[i]|shift + a * x[i] + a * x[i]*x[i], yerror[i]);

}

}
"""

NP=2;
data = {'N': len(x), 'x': x, 'y': y,'yerror':yerror,'NP':NP}
sm = pystan.StanModel(model_code=model)
fit = sm.sampling(data=data, iter=2000, chains=4, warmup=400, thin=3, seed=101)
print(fit)
``````

The output is:

``````INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_1cfdf39512ef7dc674a1650663817e89 NOW.
WARNING:pystan:Rhat above 1.1 or below 0.9 indicates that the chains very likely have not mixed
WARNING:pystan:640 of 2136 iterations saturated the maximum tree depth of 10 (30 %)
WARNING:pystan:Run again with max_treedepth larger than 10 to avoid saturation
WARNING:pystan:Chain 1: E-BFMI = 0.0067
WARNING:pystan:Chain 3: E-BFMI = 0.00587
WARNING:pystan:Chain 4: E-BFMI = 0.0581
WARNING:pystan:E-BFMI below 0.2 indicates you may need to reparameterize your model
Inference for Stan model: anon_model_1cfdf39512ef7dc674a1650663817e89.
4 chains, each with iter=2000; warmup=400; thin=3;
post-warmup draws per chain=534, total post-warmup draws=2136.

mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
shift   0.49    0.13   0.41  -0.13   0.38   0.42   0.46   2.05     10   1.49
a    1.47    0.07   0.24   0.21   1.51   1.52   1.53   1.65     13   1.29
a    -0.2  6.0e-3   0.02  -0.21   -0.2   -0.2   -0.2  -0.07     15   1.23
lp__   -1490  1745.7 8441.7 -4.4e4  -55.9 -54.38 -53.61 -53.04     23   1.13

Samples were drawn using NUTS at Wed Nov 20 17:19:09 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).
``````

I am confused as I don;t know when to trust pystan now. I don’t know why such a small modelling change causes so much difference in the fitting. Again I am most concerned not with the different results but because Rhat does not converge to 1.0 in the latter case.