Hello,

I have been trying to use Stan (via the Julia interface) for my thesis work but I keep running into the problem that the chains are all single valued. I think it has something to do with the large range of values with my inputs but I can’t be sure since this is the first real problem I have tried to to solve using Stan. I have listed my model and run summary below.

Any help would be appreciated.

## Stan Summary:

```
Inference for Stan model: model
4 chains: each with iter=(1000,1000,1000,1000); warmup=(0,0,0,0); thin=(1,1,1,1); 4000 iterations saved.
Warmup took (565, 392, 432, 356) seconds, 29 minutes total
Sampling took (70, 184, 162, 195) seconds, 10 minutes total
Mean MCSE StdDev 5% 50% 95% N_Eff N_Eff/s R_hat
lp__ -3.0e+15 3.1e+13 1.9e+15 -6.3e+15 -1.8e+15 -1.7e+15 4000 6.5e+00 7.7e+05
accept_stat__ 5.0e-01 1.4e-04 8.8e-03 5.0e-01 5.0e-01 5.0e-01 4000 6.5e+00 1.0e+00
stepsize__ 1.8e-08 5.8e-09 8.2e-09 8.2e-09 2.0e-08 3.1e-08 2.0 3.3e-03 3.1e+14
treedepth__ 7.5e-04 5.6e-04 3.5e-02 0.0e+00 0.0e+00 0.0e+00 4000 6.5e+00 1.0e+00
n_leapfrog__ 1.0e+00 7.1e-04 4.5e-02 1.0e+00 1.0e+00 1.0e+00 4000 6.5e+00 1.0e+00
divergent__ 1.0e+00 2.5e-04 1.6e-02 1.0e+00 1.0e+00 1.0e+00 4000 6.5e+00 1.0e+00
energy__ 3.0e+15 3.1e+13 1.9e+15 1.7e+15 2.0e+15 6.3e+15 4000 6.5e+00 1.5e+05
F[1] -5.5e-02 5.0e-01 7.1e-01 -1.2e+00 3.1e-01 6.7e-01 2.0 3.3e-03 4.8e+08
F[2] 8.3e-02 7.3e-01 1.0e+00 -1.4e+00 3.0e-01 1.4e+00 2.0 3.3e-03 4.7e+14
F[3] -1.6e-02 4.2e-01 6.0e-01 -7.7e-01 4.8e-01 6.5e-01 2.0 3.3e-03 6.1e+14
F[4] -2.0e-01 7.4e-01 1.0e+00 -1.3e+00 2.7e-01 1.3e+00 2.0 3.3e-03 1.2e+15
F[5] -7.2e-02 3.6e-01 5.1e-01 -9.2e-01 1.5e-01 4.2e-01 2.0 3.3e-03 2.2e+06
F[6] 2.5e-01 1.5e-01 2.2e-01 -1.7e-02 3.0e-01 5.7e-01 2.0 3.3e-03 3.3e+05
F[7] -2.5e-01 6.2e-01 8.8e-01 -1.6e+00 4.0e-01 6.8e-01 2.0 3.3e-03 1.5e+06
F[8] -8.6e-02 3.8e-01 5.4e-01 -5.6e-01 9.2e-03 7.6e-01 2.0 3.3e-03 4.7e+06
F[9] -8.2e-02 4.0e-01 5.6e-01 -7.4e-01 1.7e-01 7.1e-01 2.0 3.3e-03 6.4e+14
```

## Stan Model

```
data {
int<lower=1> N; // Number of Data Points = 8883
int<lower=1> K; // Number of Parameters = 500 but could be up to 4000
real rho; // range [0.01,0.3]
real sig_F; // 1e16
vector[N] d; // extrema(d) = (-2.2e13, 9.7e14)
vector[N] err; // extrema(err) = (6.2e12, 2.8e13)
vector[K] x; // range [0, 1]
vector[K] y; // range [0, 1]
vector[K] z; // range [0, 1]
matrix[N,K] W; // extrema(W[W .> 0]) = (4.17e-23,9e-2), mean(W[W .> 0]) = 2.8e-4
}
transformed data {
vector[K] mu;
matrix[N,K] A;
real matrix_min;
cov_matrix[K] Sigma_F;
// Scale Transfer matrix to avoid initialization error
matrix_min = max(W);
for (i in 1:K){
for (j in 1:N){
if ( (W[j,i] != 0.0) && (W[j,i] < matrix_min)){
matrix_min = W[j,i];
}
}
}
for (i in 1:K){
A[:,i] = (1/matrix_min)*W[:,i];
}
mu = rep_vector(0.0,K);
// Create Covariance Matrix
for (i in 1:K){
for (j in i:K){
Sigma_F[j,i] = ((sig_F*matrix_min)^2)*exp(-0.5*((x[i] - x[j])^2 + (y[i] - y[j])^2 + (z[i] - z[j])^2)/rho^2);
Sigma_F[i,j] = Sigma_F[j,i];
}
}
}
parameters {
//real sig_F;
//real rho; would also like to infer these
vector[K] F;
}
model {
F ~ multi_normal(mu, Sigma_F);
d ~ normal(A*F,err);
}
generated quantities {
vector[K] X;
X = F/matrix_min;
}
```