Custom pearson type 3 distribution

As a student, I often apply the Pearson type III distribution, a variant of the gamma distribution, in my academic field:

Unfortunately, the brms package does not include this distribution in the brms families. I followed the tutorial to write the code for the Pearson type III distribution:

library(brms)
library(lmomco)

log_lik_pearson3 <- function(i, prep) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  sigma <- brms::get_dpar(prep, "sigma", i = i)
  skew <- brms::get_dpar(prep, "skew", i = i)
  y <- prep$data$Y[i]
  pdfpe3(y, vec2par(c(mu, sigma, skew), type = 'pe3'))
}

posterior_predict_pearson3 <- function(i, prep, ...) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  sigma <- brms::get_dpar(prep, "sigma", i = i)
  skew <- brms::get_dpar(prep, "skew", i = i)
  quape3(runif(prep$nprep), vec2par(c(mu, sigma, skew), type = 'pe3'))
}

posterior_epred_pearson3 <- function(prep) {
  mu <- brms::get_dpar(prep, "mu")
  return(mu)
}

pearson3 <- function(link = "identity", link_sigma = "identity", link_skew = 'identity'){
  custom_family(
    "pearson3",
    dpars = c("mu", "sigma", 'skew'),
    links = c(link, link_sigma, link_skew),
    lb = c(NA, 0, 0),
    ub = c(NA, NA, NA),
    type = "real",
    # log_lik = log_lik_pearson3,
    # posterior_predict = posterior_predict_pearson3,
    # posterior_epred = posterior_epred_pearson3
  )
}

stan_pearson3 <- "
  real pearson3_lpdf(real y, real mu, real sigma, real skew) {
    real shape = (2 / skew)^2;
    real location = sqrt((sigma^2 / shape));
    real scale = mu - shape * location;
    real x = (y - scale) / location;

    return - log(location) - 
           log(tgamma(shape)) + (shape - 1) * log(x) -
           x;
  }
"
stanvars = stanvar(scode = stan_pearson3, block = "functions")

fit1 <- brm(formula = mpg ~ wt,
            family = pearson3, stanvars = stanvars, data = mtcars,
            control = list(adapt_delta = .90), iter = 2000)

summary(fit1); plot(fit1)

but the results were too inefficient:

Warning messages:
1: There were 1796 divergent transitions after warmup. See
Runtime warnings and convergence problems
to find out why this is a problem and how to eliminate them.
2: Examine the pairs() plot to diagnose sampling problems
3: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
Runtime warnings and convergence problems

Moreover, I couldn’t find errors or make improvements within my ability. I would be very grateful if you had time to check it out!

Thanks!