I am tyring to create a custom family in brms using an lpdf that works fine in rstan. When I try to run the code below, I get an error that says “could not find function vreal”. I’m guessing that the problem is in the custom_family function below or in the specification in the brm() model, but I cannot find the solution to it anywhere else online or on this forum.

```
paretocounts <- custom_family(
"paretocounts", dpars = c("mu"),
links = c("identity"),
lb = -Inf, ub = Inf,
type = "real", vars = c("vint1[counts]",
"vreal1[xmin]",
"vreal2[xmax]")
)
stan_funs <- "
real paretocounts_lpdf(real x, real mu, real xmin, real xmax, real counts){
if(mu != -1)
return(counts*(log((mu+1) / ( xmax^(mu+1) - xmin^(mu+1))) + mu*log(x)));
else
return(counts*(log(log(xmin) - log(xmax)) + mu*log(x)));
}
"
stanvars <- stanvar(scode = stan_funs, block = "functions")
#### load attached data called dw_short.csv ####
fit2 <- brm(
dw | vint(counts) | vreal(xmin, xmax) ~ mat_s + (1|sample),
data = dw_short,
family = paretocounts,
stanvars = stanvars
)
```

Data attached

Please also provide the following information in addition to your question:

- Operating System: Windows 10 Enterprise
- brms Version: 2.19.0

dw_short.csv (494.4 KB)

just don’t use vint and vreal in the vars argument of custom_family(). brms handles that itself.

Thanks for the quick reply. How should the addition terms be added to the brm formula? Currently I have them as follows:

fit2 ← brm(

dw | vint(counts) | vreal(xmin, xmax) ~ mat_s + (1|sample),

data = dw_short,

family = paretocounts,

stanvars = stanvars

)

That didn’t work to begin with. I’ve tried it by removing the addition terms and the error message appears to be looking for the addition terms somewhere:

Error in stanc(file = file, model_code = model_code, model_name = model_name, :

0

Semantic error in ‘string’, line 58, column 16 to column 47:

Ill-typed arguments supplied to function ‘paretocounts_lpdf’. Available signatures:

(real, real, real, real, real) => real

Instead supplied arguments of incompatible type: real, real.

they should be there but don’t use the vint part of the vars argument. also combine multiple addition terms with a +. please see the brms_customfamilies vignette for more details

Thank you. It’s getting closer but still gives an error and I couldn’t find the solution in the vignette. I think I’m confused about how to enter *multiple* addition terms in the formula. The current error is this:

Error: The following addition terms are invalid:

‘counts’, ‘xmin’, ‘xmax’

Based on this code:

paretocounts ← custom_family(

“paretocounts”, dpars = c(“mu”),

links = c(“identity”),

lb = -Inf, ub = Inf,

type = “real”, vars = c(“counts”,

“xmin”,

“xmax”))

stan_funs ← "

real paretocounts_lpdf(real Y, real mu, real xmin, real xmax, real counts){

if(mu != -1)

return(counts*(log((mu+1) / ( xmax^(mu+1) - xmin^(mu+1))) + mu*log(Y)));*

else

return(counts(log(log(xmin) - log(xmax)) + mu*log(Y)));

}

"

stanvars ← stanvar(scode = stan_funs, block = “functions”)

dw_short = dw %>% filter(sample <= 4)

fit2 ← brm(

dw | (counts + xmin + xmax) ~ mat_s + (1|sample),

data = dw_short,

family = paretocounts,

stanvars = stanvars

)

But that gives an error that:

Thank you. If there are two vreals, how can we add them both?

This doesn’t seem to work:

y|vint(name1) + vreal(name2) + vreal(name3) ~ …

I get the error that addition terms can only occur once.

please read the docs more carefully. you are asking many questions that can be inferred from them. you do vreal(x1,x2)

Thanks. I’ll keep working on it through the docs. I appreciate the help. It is much closer than when I started.

Solved.

```
rparetocounts <- function(n = 300, mu = -1.2, vreal2 = 1, vreal3 = 1000) {
samples <- numeric(n)
{
if(vreal2 <= 0 | vreal2 >= vreal3) stop("Parameters out of bounds in rPLB")
u <- stats::runif(n)
if(mu != -1)
{ y <- ( u*vreal3^(mu+1) + (1-u) * vreal2^(mu+1) ) ^ (1/(mu+1))
} else
{ y <- vreal3^u * vreal2^(1-u)
}
return(y)
}
return(samples)
}
dparetocounts <- function(x, mu, vreal2, vreal3) {
if (vreal2 <= 0 || vreal2 >= vreal3)
stop("Parameters out of bounds in dPLB")
if (x < vreal2 || x > vreal3)
return(0)
if (mu != -1) {
density <- (mu + 1) * (x^(mu+1)) / (vreal3^(mu+1) - vreal2^(mu+1))
} else {
density <- x^(-2) / (vreal2 * log(vreal3/vreal2))
}
density
}
log_lik_paretocounts <- function(i, prep) {
mu <- brms::get_dpar(prep, "mu", i = i)
y <- prep$data$Y[i]
dparetocounts(x, mu, vreal2, vreal3)
}
posterior_predict_paretocounts <- function(i, prep, ...) {
mu <- brms::get_dpar(prep, "mu", i = i)
vreal2 = prep$data$vreal2[i]
vreal3 = prep$data$vreal3[i]
rparetocounts(prep$ndraws, mu, vreal2, vreal3)
}
posterior_epred_paretocounts <- function(prep) {
mu <- prep$dpars$mu
return(mu)
}
stan_funs <- "
real paretocounts_lpdf(real Y, real mu, real vreal1, real vreal2, real vreal3){
if(mu != -1)
return(vreal1*(log((mu+1) / ( vreal3^(mu+1) - vreal2^(mu+1))) + mu*log(Y)));
else
return(vreal1*(log(log(vreal2) - log(vreal3)) + mu*log(Y)));
}
"
stanvars <- brms::stanvar(scode = stan_funs, block = "functions")
paretocounts <- function(){brms::custom_family(
"paretocounts",
dpars = c("mu"),
links = c("identity"),
lb = -Inf,
ub = Inf,
type = "real",
vars = c("vreal1[n]",
"vreal2[n]",
"vreal3[n]"),
posterior_predict = posterior_predict_paretocounts,
posterior_epred = posterior_epred_paretocounts,
log_lik = log_lik_paretocounts)
}
```