Beta_proportion parameter recovery

Modeling
,
1

Hey everyone,

I’ve tried to start working with beta_proportion sampling for a constrained parameter [0,1]. I started with a toy model where I hierarchically sampled proportions for a Bernoulli process and then tried to recover the proportions. However, when I sampled values from a beta distribution with a rather small sd (high kappa) it led to difficulties in precisely recovering it (big sd).

variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ -121152.13 -121152.00 24.91 25.20 -121193.00 -121113.00 1.00 1223 2099
mu 0.30 0.30 0.00 0.00 0.30 0.30 1.00 3795 3046
kappa 86.46 86.34 5.48 5.37 77.82 95.69 1.00 2318 2669

I attach the stan and R code snippets.

data {

  int<lower=0> N;

  

  int<lower=0> trials;

  

  array[N, trials] int<lower=0> choice;

}

parameters {

  real<lower=0, upper=1> mu;

  

  real<lower=0> kappa;

  

  vector<lower=0, upper=1>[N] alpha;

}

model {

  mu ~ beta(2, 2);

  

  kappa ~ uniform(0, 200);

  

  for (i in 1 : N) {

    alpha[i] ~ beta_proportion(mu, kappa);

    for (t in 1 : trials) {

      choice[i, t] ~ bernoulli(alpha[i]);

    }

  }

}

my_compiledmodel <- cmdstan_model(path,compile=T)

# Prepare data for Stan model
# Number of individuals
N <- 1000  

# Number of observations per individual
trials <- 200
# Desired mean and standard deviation
mean = 0.3
sd = 0.05
kappa = (mean*(1-mean)/sd^2)-1

# Calculate alpha and beta from mean and sd
alpha = mean * ((mean * (1 - mean) / sd^2) - 1)
beta = (1 - mean) * ((mean * (1 - mean) / sd^2) - 1)

# Generate random numbers from the beta distribution
random_numbers = rbeta(N, alpha, beta)

# Generate data
y <- matrix(nrow = N, ncol = trials)
for (i in 1:N) {
  y[i, ] <- rbinom(trials, 1, random_numbers[i])
}

df=as.data.frame(y)
# Prepare a list of data for Stan
stan_data <- list(N = N, trials = trials, choice = y)


# Fit the model
fit<- my_compiledmodel$sample(
  data            = stan_data, 
  iter_sampling   = 1000,
  iter_warmup     = 1000,
  chains          = 4,
  parallel_chains = 4,
  refresh         = 20)

# Extract the results
fit

Thanks a lot!