Custom rng function beta distribution

Modeling
, ,
1

Hi,

How can I specify the rng function for the beta distribution? I know that it is implemented in Stan, but I would like to know how to specify it manually.

thank you in advance.

2

The nice thing about Stan is that it’s all available to see math/stan/math/prim/prob/beta_rng.hpp at develop · stan-dev/math · GitHub.

Using chatgpt I converted that to

log_sum_exp <- function(a, b) {
  max_val <- max(a, b)
  log(exp(a - max_val) + exp(b - max_val)) + max_val
}

beta_rng <- function(alpha, beta) {
 
  N <- max(length(alpha), length(beta))
  output <- numeric(N)
  
  for (n in 1:N) {
    if (alpha[n] > 1.0 && beta[n] > 1.0) {
      a <- rgamma(1, alpha[n], 1.0)
      b <- rgamma(1, beta[n], 1.0)
      output[n] <- a / (a + b)
    } else {
      log_a <- log(runif(1)) / alpha[n] + log(rgamma(1, alpha[n] + 1, 1.0))
      log_b <- log(runif(1)) / beta[n] + log(rgamma(1, beta[n] + 1, 1.0))
      log_sum <- log_sum_exp(log_a, log_b)
      output[n] <- exp(log_a - log_sum)
    }
  }
  
  return(output)
}

Let’s test and see if we can return the correct values

data {
  int<lower=0> N;
  vector[N] y;
}
parameters {
  real<lower=0> alpha;
  real<lower=0> beta;
}
model {
  y ~ beta(alpha, beta);
}

And call it in R

x <- c()
for(i in 1:500) 
  x[i] <- beta_rng(0.4, 2.4)

library(cmdstanr)
mod <- cmdstan_model("beta.stan")

mod_out <- mod$sample(
  data = list(N = 500,
              y = x),
  parallel_chains = 4
)
mod_out

With output

 variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
    lp__  660.19 660.51 1.02 0.73 658.15 661.16 1.00     1934     2275
    alpha   0.40   0.40 0.02 0.02   0.36   0.43 1.00     1727     2049
    beta    2.46   2.46 0.19 0.19   2.16   2.78 1.00     1837     2371
4 Likes