@Bob_Carpenter thank you so much for your comment. I read the Stan manual, but I could not find it helpful to the problem above, because I want to see a sample code that demonstrates the use of Bayesian inference (your explanation is much better). I tried modifying the original code above, and adding the main R files that executes the â€śmodelStringâ€ť below, but the code failed to run (it keeps giving me the error message: PARSER EXPECTED: whitespace to end of file. FOUND AT line 35: }, but I checked my code and saw no whitespace at the end of line 35). Also, I am not sure if the code does what i want it do though, so could you please help review the code below to see if my coding setup is correct? I really appreciate your help!

```
library("truncnorm")
modelString= "data {
# // Define log probability density function
# real myNormal_lpdf(real x, real A%*%yT, real sigma2*rep(1, rtruncnorm(ncol, 0, 1))) {
# return -log(2 * pi()) / 2 - log(sigma2*rep(1, rtruncnorm(ncol, 0, 1)))
# - square(A%*%yT - x) / (2 * (sigma2*rep(1, rtruncnorm(ncol, 0, 1)))^2);
# }
# }
# data {
# int N;
# real yh[N];
}
# parameters {
# real<lower = 0> sigma2;
# real<lower=0> sigma;
real<lower = 0> x;
# }
# model {
# // Priors
# yT ~ normal(yh, sigma*rep(1, rtruncnorm(ncol, 0, 1)));
# // Likelihood
# for(n in 1:N) {
# target += myNormal_lpdf(x| A%*%yT, sigma2*rep(1, rtruncnorm(ncol, 0, 1)));
# }
# }
# generated quantities{
# real yh_ppc;
# {
# real x;
# x = uniform_rng(0, 1);
# // Phi is the probit function in Stan, the CDF of the standardised Normal N(0,1)
# // inv_Phi is the quantile function of the standardised Normal N(0,1)
# yh_ppc = mu + sigma * inv_Phi(x);}
}"
library(rstan)
library(MASS)
library(truncnorm)
set.seed(1)
nrow=16;
ncol = 24;
sigma=0.2;
sigma2=0.3;
lambda = matrix(sample.int(100, size = ncol*nrow, replace = T),nrow,ncol);
lambda = lambda - diag(lambda)*diag(x=1, nrow, ncol);
yh = rpois(ncol,lambda) + rtruncnorm(ncol,0,1,mean = 0, sd = 1);
A = matrix(c(1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_A^1 - Done
0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_A^2 - Done
0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_B^1 - Done
1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_B^2 - Done
0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_C^1 - Done
0,1,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0, #x_C^2 - Done
0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0, #x_D^1 - Done
0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_D^2 - Done
0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0, #x_E^1 - Done
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,1,0, #x_E^2 - Done
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0, #x_F^1 - Done
0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,1,0,0, #x_F^2 - Done
0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0, #x_G^2 - Done
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1, #x_G^1 - Done
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0, #x_H^1 - Done
0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,1), #x_H^2 - Done
nrow, ncol, byrow= TRUE);
truehist(x, col="#B2001D")
lines(density(x), col="skyblue", lwd=2)
summary(yT)
ret_sm <- stan_model(model_code = modelString) # Compile Stan code
fit <- sampling(ret_sm, warmup=100, iter=600, seed=1,
data=list(yh, N=length(yh)))
stan_trace(fit, inc_warmup = TRUE)
stan_hist(fit)
print(fit, probs=c(0.025, 0.5, 0.975))
summary(extract(fit, "yh_ppc")[["yh_ppc"]])
plot(ecdf(x), main="Posterior predictive check")
lines(ecdf(extract(fit, "yh_ppc")[["yh_ppc"]]), col="#B2001D")
```