# Custom rng function beta distribution

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.

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