I’m trying to implement P-splines with PC prior in stan and results are less than stellar.

In particular:

```
The E-BFMI, 0.00, is below the nominal threshold of 0.30 which suggests that HMC may have trouble exploring the target distribution.
```

Also a lot of parameters have Rhat>1.05 and fewer than 0.001 effective draws per transition.

So I’m obviously doing something bad.

```
data {
int<lower=0> N; // number of data items
int<lower=0> K; // number of splines
int<lower=0> L; // number of simulated samples
matrix[N, K] B; // predictor matrix
matrix[K, K] R; // structure matrix
vector[N] y; // outcome vector
real U; // upper bound on marginal std
real alpha; //probability of marginal std exceeding U
}
transformed data {
real theta = -log(alpha)/U;
vector[L] ones_L = rep_vector(1, L);
vector[K] zeros_K = zeros_vector(K);
matrix[L,K] sub_B = block(B,1,1,L,K);
}
parameters {
vector[K] beta; // coefficients
real<lower=0> sigma; // std
real<lower=0> tau_beta; // beta precision
}
model {
y ~ normal(B * beta, sigma); // likelihood
beta ~ multi_normal_prec(zeros_K,tau_beta*R);
tau_beta ~ gumbel(0.5,theta); //pc prior on precision of beta
sigma ~ exponential(.1);
}
generated quantities {
array[L] real y_pred; //predictions
y_pred = normal_rng(sub_B*beta,ones_L*sigma);
}
```

I specified parameters as:

```
nm1= BN.shape[1]-1
D=np.concatenate((np.zeros((nm1,1)),-np.identity(nm1)),axis=1)+np.concatenate((np.identity(nm1),np.zeros((nm1,1))),axis=1)
R = D.T@D
data_current = dict(N=BN.shape[0],
K=BN.shape[1],
L=M,
B=BN,
R=R,
y=melted_traj.current.values[0:N*M],
U=5,
alpha=0.05)
```

BN is a stacked matrix of predictions for multiple signals (10).

I’m using PC-priors description from: [1511.05748] Penalized complexity priors for degrees of freedom in Bayesian P-splines

I’d be greatful for suggestion what I am doing wrong.