Using the beta_proportion distribution of Stan

stan_beta_proportion.Rmd

---
title: "R Notebook"
output: html_notebook
---

```{r}
library(rstan)
options(mc.cores = parallel::detectCores())
library(bayesplot)
```

Data generation

```{r data}
inv_logit <- function(x) {
  1 / (1 + exp(-x))
}

set.seed(123)
N <- 100
X <- runif(N, -1, 1)
mu <- inv_logit(-1 + 2 * X)
kappa <- 2
alpha <- mu * kappa
beta <- (1 - mu) * kappa
Y <- rbeta(N, alpha, beta)
ggplot(data.frame(X = X, Y = Y)) +
  geom_point(aes(X, Y))
```

Stan model

```{stan output.var=beta_prop}
data {
  int<lower = 0> N;
  vector[N] X;
  vector<lower = 0, upper = 1>[N] Y;
}
parameters {
  vector[2] beta;
  real<lower = 0> kappa;
}
transformed parameters {
  vector[N] mu = inv_logit(beta[1] + beta[2] * X);
}
model {
  Y ~ beta_proportion(mu, kappa);
}
generated quantities {
  vector[N] yrep;

  for (n in 1:N)
    yrep[n] = beta_proportion_rng(mu[n], kappa);
}
```

Fit

```{r fit}
data <- list(N = N, X = X, Y = Y)
fit <- sampling(beta_prop, data = data, iter = 2000)
print(fit, pars = c("beta", "kappa"))
```

Posterior predictive check
```{r ppc}
yrep <- extract(fit, pars = "yrep")[[1]]
ppc_dens_overlay(Y, yrep[sample(nrow(yrep), 100), ])
```