# Brms custom family continuous bernoulli

Hi all -

In the blog post below I discuss the continuous Bernoulli distribution and how to fit it in brms as a custom family:

Inspired by code from @spinkey. To see just the R code to fit the model, I copied and pasted it below:

library(ordbetareg)
library(brms)

N <- 500
X <- runif(n=N)
Y <- rordbeta(n=N, mu = plogis(-2 + 2.5*X),cutpoints=c(-3,3))

hist(Y)

# define custom family

c_bernoulli <- custom_family("c_bernoulli",
dpars="mu",
lb=0,ub=1,
type="real")

# define log density function
# some code from spinkey https://discourse.mc-stan.org/t/continuous-bernoulli/26886

stan_funs <- "

//normalization constant
real c_norm(real mu) {

if(mu==0.5) {

return log(2);

} else {

real const = (log(2 - 2*mu) - log(2*mu))/(2 * (1 - 2*mu));
return(log(const));

}

}

// log PDF for continuous Bernoulli
real c_bernoulli_lpdf(real y, real mu) {

// unnormalized density

real lp = y * log(mu) + (1 - y) * log1m(mu);

// normalized density

lp += c_norm(mu);

return lp;

}"

stanvars <- stanvar(scode = stan_funs, block = "functions")

# posterior predictions

posterior_predict_c_bernoulli <- function(i, prep, ...) {

# need inverse CDF function for continuous bernoulli

inv_cdf <- function(u=NULL, mu=NULL) {

if(mu==0.5) {

out <- u

} else {

out <- (log(u * (2 * mu - 1) + 1 - mu) - log(1 - mu))/(log(mu) - (log(1-mu)))

}

return(out)

}

mu <- brms::get_dpar(prep, "mu", i = i)
u <- runif(n=length(mu))

inv_cdf(u,mu)

}

posterior_epred_c_bernoulli <- function(prep) {

# expected value
mu <- brms::get_dpar(prep, "mu")

if(mu==0.5) {

return(0.5)

} else {

(mu / (2 * mu - 1)) + (1 / (2*atanh(1 - 2*mu)))

}

}

fit_c_bern <- brm(
Y ~ X, data = tibble(Y=Y, X=X),
family = c_bernoulli, stanvars = stanvars, backend="cmdstanr"
)

Would you be up to open a PR at GitHub - sims1253/bayesfam to add it? Looks cool :)

At some point I’ll get around to adding it to ordbetareg, so can consider at that point adding it to your github.