Hi everyone, i just came across this post! Not sure if it is appropriate to ask it here but ill ask anyway. I am currently a university student working on one of Pierre’s Jacob’s paper “Better together? Statistical learning in models made of modules” for my thesis. I am attempting to code the cut posterior of the biased data on Stan. So far i have this and it doesnt seem to work! (My prof die die wants me to use it on Stan)

```
library(StanHeaders)
library(Rcpp)
library(ggplot2)
library(rstan)
library(withr)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
N1 <-100
N2 <-1000
Y1 <- rnorm(N1, mean = 0, sd = 1)
Y2 <- rnorm(N2, mean = 1, sd = 1)
lambda1 <- 1
lambda2 <- 100
model_data1 <- list(N1=N1,Y1 = Y1,lambda1 = lambda1, lambda2 = lambda2)
fit2 <-stan(file = 'model1.stan', data = model_data1)
```

```
//model1.stan
data
{ int N1;
vector[N1] Y1;
real lambda1;
}
parameters {
real theta1;
}
model {
for(i in 1:N1)
Y1[i] ~ normal(theta1, 1); theta1 ~ normal(0,lambda1^-1);
}
```

```
my_samples2 <- extract(fit2)
theta1samples <- my_samples2$theta1[2001:2100]
model_data4<-list(N1=N1,N2=N2,Y2 = Y2,theta1samples=theta1samples,lambda2 = lambda2)
fit5 <- stan(file = 'cut1ci.stan', data = model_data4)
my_samples4 <- extract(fit5)
```

```
//cut1ci.stan
data {
int N1;
intN2;
vector[N2] Y2;
real lambda2;
vector[N1] theta1samples;
}
parameters {
real theta2;
}
model {
for(i in 1:N1)
Y2[i] ~ normal(theta1samples[i] + theta2, 1);
theta2 ~ normal(0, sqrt(inv(lambda2)));
}
```

Then, if both draws are correct, how do i put both draws together? Thank you for your help!

**ADMIN EDIT:** Moved to a new thread and added code formatting by @martinmodrak