[edit: escaped Stan program]

Hello, here I am trying to fit a dual-rate motor adaptation model using rstan. The model has 6 parameters, retention rates: Af, As. Learning rates: Bf, Bs. And state noise and execution noise standard deviation.

It also has 2 hidden states - a fast state and a slow state, the constraint is such that As > Af and Bf > Bs, and As, Af, Bs, Bf belongs to (0,1). The noise are believed to be gaussian center at 0, with std >0.

y is the observed data, and p-y is the error term associated with learning rate, nofeedback is just an indictor for whether the error is present.

Here is my model:

```
data {
int<lower=1> NSubjects;
int<lower=1> TrialEnd[NSubjects];
int<lower=1> max_length;
matrix[max_length, NSubjects] y;
matrix[max_length, NSubjects] p;
matrix[max_length, NSubjects] nofeedback;
}
parameters {
vector<lower=0, upper=1>[NSubjects] Af;
vector<lower=0, upper=1>[NSubjects] Bs;
vector<lower=0, upper=1>[NSubjects] delta_A; // Ensuring a minimum difference
vector<lower=0, upper=1>[NSubjects] delta_B; // Ensuring a minimum difference
vector<lower=0>[NSubjects] execute_noise;
vector<lower=0>[NSubjects] state_noise;
}
transformed parameters {
vector<lower=0, upper=1>[NSubjects] As;
vector<lower=0, upper=1>[NSubjects] Bf;
for (i in 1:NSubjects) {
As[i] = Af[i] + delta_A[i]; // As is now guaranteed to be larger than Af
Bf[i] = Bs[i] + delta_B[i]; // Bf is now guaranteed to be larger than Bs
}
}
model {
// prior
Af ~ normal(0.7, 0.1);
Bs ~ normal(0.1, 0.1);
delta_A ~ normal(0.2, 0.05);
delta_B ~ normal(0.2, 0.05);
execute_noise ~ cauchy(0, 2);
state_noise ~ cauchy(0, 2);
for (j in 1:NSubjects) {
vector[max_length] fast;
vector[max_length] slow;
fast[1] = 0;
slow[1] = 0;
for (n in 2:TrialEnd[j]) {
fast[n] = Af[j] * fast[n-1] + Bf[j] * (p[n-1,j] - y[n-1,j]) * (1 - nofeedback[n-1,j]);
slow[n] = As[j] * slow[n-1] + Bs[j] * (p[n-1,j] - y[n-1,j]) * (1 - nofeedback[n-1,j]);
y[n,j] ~ normal(fast[n] + slow[n], sqrt(execute_noise[j]^2 + state_noise[j]^2));
}
}
}
```

The problem is parameters dont converge well (no rhat <= 1.05) also I see a high correlation between Af and As, Bf and Bs(could be caused by the way I parameterize the constraint). Is there any way to make the model more robust? Thank you !!