Hello folks,
I am trying to implement the Generalized asymmetric Laplace (GAL) as a likelihood in Stan using brms.
Link to paper that introduce and described the GAL distribution here
I have the following implementation. It works for some quartile but tends to generate a lot of warnings for the median.
1: There were 838 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
The code below is what I have to help reproduce the issue. You will note that for 0.25 or 0.75, the code runs flawlessly. However, for quantile tau0 = 0.5 (or around 0.5), It generates lost of warnings.
Does any one have suggestions or hint on what’s causing this?
#--- This code will illustrate some of the issues with implementation
#--- GAL in STAN
library(brms)
#-------------------------------------------------------------------------------
#--- Adhoc functions
GamF <- function(gam, p0){
(2*pnorm(-abs(gam))*exp(.5*gam^2) - p0)^2
}
#optimize(GamF,interval = c(-30,30), p0=1 - 0.05)
GamBnd <- function(p0){
Re1 <- optimize(GamF, interval = c(-30, 30), p0 = 1-p0)
Re2 <- optimize(GamF, interval = c(-30, 30), p0 = p0)
c( -abs(Re1$minimum), abs(Re2$minimum), Re1$objective,Re2$objective)
}
#-------------------------------------------------------------------------------
#-------- Sample a data
#-------------------------------------------------------------------------------
set.seed(0)
N = 500; X = cbind(rbinom(n=N, size = 1, prob = 0.65), runif(n=N, -1,1)); beta = c(-1,2); Y = 1 + X%*%beta + rnorm(n=N, sd=1)
Data <- data.frame(Y =Y, X1 = X[,1], X2 = X[,2])
#--- Function for the GAL
tau0 = 0.5 # here you can change the quantile level between (0, 1)
Bd = GamBnd(tau0)[1:2]
{
GAL2 <- custom_family(
"GAL2", dpars = c("mu", "sigma","gam", "tau"), links=c("identity","log","identity","identity"),
lb=c(NA,0, Bd[1]*0.95,0), ub=c(NA,NA, Bd[2]*0.95,1), type="real") #, vars = "vint1[n]"
stan_funs2 <- "
/*
A = -est*p_neg + .5*pow(gam, 2)*pow(p_neg/p_pos, 2) + log(Phi_approx(a2-a3)) + log1m_exp(fabs(log(Phi_approx(a2-a3)) - log(Phi_approx(a2))));
gam = (gamU - gamL) * ligam + gamL;
real gam = (gamU - gamL) * ligam + gamL;
real GAL2_lpdf(real y, real mu, real sigma, real ligam, real tau, real gamL, real gamU){
real GAL2_rng(real mu, real sigma, real ligam, real tau, real gamL, real gamU){
*/
/* helper function for asym_laplace_lpdf
* Args:
* y: the response value
* tau: quantile parameter in (0, 1)
*/
real rho_quantile(real y, real tau) {
if (y < 0) {
return y * (tau - 1);
} else {
return y * tau;
}
}
/* asymmetric laplace log-PDF for a single response
* Args:
* y: the response value
* mu: location parameter
* sigma: positive scale parameter
* tau: quantile parameter in (0, 1)
* Returns:
* a scalar to be added to the log posterior
*/
real asym_laplace_lpdf(real y, real mu, real sigma, real tau) {
return log(tau * (1 - tau)) -
log(sigma) -
rho_quantile((y - mu) / sigma, tau);
}
/* asymmetric laplace log-CDF for a single quantile
* Args:
* y: a quantile
* mu: location parameter
* sigma: positive scale parameter
* tau: quantile parameter in (0, 1)
* Returns:
* a scalar to be added to the log posterior
*/
real asym_laplace_lcdf(real y, real mu, real sigma, real tau) {
if (y < mu) {
return log(tau) + (1 - tau) * (y - mu) / sigma;
} else {
return log1m((1 - tau) * exp(-tau * (y - mu) / sigma));
}
}
/* asymmetric laplace log-CCDF for a single quantile
* Args:
* y: a quantile
* mu: location parameter
* sigma: positive scale parameter
* tau: quantile parameter in (0, 1)
* Returns:
* a scalar to be added to the log posterior
*/
real asym_laplace_lccdf(real y, real mu, real sigma, real tau) {
if (y < mu) {
return log1m(tau * exp((1 - tau) * (y - mu) / sigma));
} else {
return log1m(tau) - tau * (y - mu) / sigma;
}
}
real GAL2_lpdf(real y, real mu, real sigma, real gam, real tau){
real p_pos;
real p_neg;
real a3;
real a2;
real p;
real est;
real A;
real B;
real Res = 0;
//real gam = ligam;
p = 1 * (gam < 0) + (tau - 1 * (gam < 0))/(2*Phi(-fabs(gam)) * exp(.5 * pow(gam, 2)));
p_pos = p - 1 * (gam > 0);
p_neg = p - 1* (gam < 0);
est = (y - mu) / sigma;
if(fabs(gam) > 0){
a3 = p_pos * (est / fabs(gam));
a2 = fabs(gam) * (p_neg / p_pos);
if(est/gam > 0){
A = 0.5 * pow(gam, 2) * pow(p_neg/p_pos, 2) - est * p_neg + log_diff_exp(log(Phi_approx(a2-a3)), log(Phi_approx(a2)) );
B = 0.5 * pow(gam, 2) - p_pos * est + log(Phi_approx(-fabs(gam) + a3));
Res = log(2*p*(1-p)) - log(sigma) + log_sum_exp(A, B);
}else{
Res = log(2*p*(1-p)) - log(sigma) - p_pos * est + 0.5 * pow(gam, 2) + log(Phi_approx(-fabs(gam) ));
}
}else{
Res = asym_laplace_lpdf( y | mu, sigma, tau);
}
return Res;
}
real GAL2_rng(real mu, real sigma, real ligam, real tau){
real A;
real B;
real C;
real p;
real hi;
real nui;
real mui=0;
real Up = uniform_rng(.5, 1.0);
real gam = ligam;
p = (gam < 0) + (tau - (gam < 0))/(2*Phi_approx(-fabs(gam))*exp(.5*pow(gam, 2)));
A = (1 - 2*p)/(p - pow(p,2));
B = 2/(p - pow(p,2));
C = 1/((gam > 0) - p);
hi = sigma * inv_Phi(Up);
nui = sigma * exponential_rng(1);
mui += mu + A * nui + C * fabs(gam) * hi;
return normal_rng(mui, sqrt(sigma*B*nui));
}
"
}
#-------------------------------------------------------------------------------
#----------------- Now run the code
#-------------------------------------------------------------------------------
Nchain = 2
Niter = 2000
formCase1b <- paste0("Y ~ X1 + X2")
stanvars2 <- stanvar(scode = stan_funs2, block = "functions")
fitCase1b <- brms::brm(bf(as.formula(formCase1b), tau = tau0), data = Data, family = GAL2, stanvars = stanvars2,
chains = Nchain,iter = Niter, control = list(adapt_delta = 0.99,max_treedepth=12),
cores = 2, seed = 1123) # init = 0.1,
print(fitCase1b)