My question relates to fitting custom distributions in R but I feel it has enough of a probability element to remain on CV.

I have an interesting set of data which has the following characteristics:

  • Large mass at zero
  • Sizeable mass below a threshold that fits a right-skewed parametric distribution very well
  • Small amount of mass at extreme values
  • A number of covariates that should drive the variable of interest

I was hoping to model this using a zero-inflated distribution approach, which is widely explored in the literature. Essentially, the density is:

$$f_{Y}(y)=\begin{cases} \pi \quad\quad\quad\quad\,\,\,\,\,\,,\,\,y=0 \\ (1-\pi)f_X(y),\,\,y>0 \end{cases}$$

This is easy enough to fit as is. However, I would like the mixing parameter $\pi$ to be dependent on the covariates $Z$ via a logistic regression:

$$\text{logit}(\mathbb{E}[\pi\,|\,Z])=\beta Z$$ where $\beta$ is a vector of coefficients for the covariates. Furthermore, because of the extreme-tail nature of my data, my distribution $f_{X}(y)$ fits best with an extreme-value approach:

$$f_{X}(y)=\begin{cases} f_{A}(y;a,b) \quad\,\,\,\,\,\,\,,\,y\leq \mu \\ (1-F_{A}(\mu))\cdot\text{GPD}\bigg(y;\mu,\sigma=\frac{(1-F_{A}(\mu))}{f_{A}(\mu)},\xi\bigg),\,y>\mu \end{cases}$$ where $\text{GPD}(y;\mu,\sigma,\xi)$ refers to the Generalized Pareto distribution, modelling the excess above a certain threshold $\mu$ and $f_{A}(y;a,b)$ is a given right-skewed distribution with scale and shape parameters $a$ and $b$, respectively. The above characterization ensures that the densities are continuous at $y=\mu$ (not differentiable, though) and that $f_{X}(y)$ integrates to 1.

In addition, I would ideally want the parameters of the above distributions to also depend on covariates:

$$f_{A}(y;a,b,\beta Z)$$ $$\text{GPD}\bigg(y;\mu,\sigma=\frac{(1-F_{A}(\mu))}{f_{A}(\mu)},\xi,\beta Z\bigg)$$

I realize that the above setup is quite complex but I was wondering if there is a way to derive the MLE estimates of each of the desired parameters by maximizing the likelihood function i.e. to obtain:

$$\hat{\xi}, \hat{a}, \hat{b}, \hat{\beta}$$

Is there an feasible/ideal way to go about this in R? Both in terms of my specific problem but also fitting custom distributions more generally?

  • 7
    $\begingroup$ The way you construct $f_X(y)$ by "cut and paste" means that $f_X(y)$ don't not integrate to one. You need to reformulate your model somehow to fix this, otherwise any estimates you might obtain will be meaningless. You also need to think about what parameters should depend on the covariates and what parameters should remain constant, perhaps after reparameterization in terms of means, variances and skews etc. of the different model parts. Apart from this, with good starting values you may be able to fit something like this using optim in R if you have enough data. $\endgroup$ Commented Aug 25, 2017 at 13:53
  • 1
    $\begingroup$ I agree with @Jarle Tufto, the $f_A(y)$ density-as specified in the post-is not guaranteed to sum/integrate to 1. The only way to ensure that happens is to choose $\mu$ accordingly and I don't think that is what the OP intended. $\endgroup$ Commented Aug 25, 2017 at 17:32
  • 1
    $\begingroup$ Is there a particular reason that you want a maximum-likelihood estimate rather than a Bayesian solution? $\endgroup$ Commented Aug 27, 2017 at 22:05
  • 3
    $\begingroup$ @StatsPlease -> do you want an algorithm that works for ANY generic $f_A(y)$ or that works for a specific distributional family, say for example, Gamma(a, b)? $\endgroup$ Commented Sep 1, 2017 at 13:13
  • 2
    $\begingroup$ @StatsPlease -> No continuous distribution without a point-mass at 0 will work with your zero-inflation specification. The reason is that you will introduce a latent variable for the zero-inflation and some of those latent variables will need to be 0 and others 1. The 0 values indicate an observed 0 from the underlying and the a 1 indicates a structural 0. In the distributions discussed in the comments (Gamma, Weibull, Pareto) all have $\Pr(Y=0)=0$. So you will have a perfect separation problem in the logit model fitting that you've asked for in the post. $\endgroup$ Commented Sep 2, 2017 at 12:59

2 Answers 2


This answer assumes $\mu$ is known.

One very flexible way to get MLE's in R is to use STAN via rstan. STAN has a reputation for being an MCMC tool, but it also can estimate parameters by variational inference or MAP. And you're free to not specify the priors.

In this case, what you're doing is very similar to their hurdle-model example. Here is the STAN code for that example.

data {
  int<lower=0> N;
  int<lower=0> y[N];
parameters {
  real<lower=0, upper=1> theta;
  real<lower=0> lambda;
model {
  for (n in 1:N) {
    if (y[n] == 0)
      target += bernoulli_lpmf(1 | theta);
      target += bernoulli_lpmf(0 | theta)
                  + poisson_lpmf(y[n] | lambda);

To adapt this for your own use, you could:

  • Replace poisson_lpmf with the log-density for your $f_A$.
  • Add a third case to the if-else so that it checks for exceeding $\mu$, not just 0. As the meat of that third case, use the log pmf for your extreme value distribution of choice.
  • Replace bernoulli_lpmf with categorical_lpmf and make the mixture probability parameter into a vector.
  • To incorporate covariates, you can add regression parameters, and make all your other parameters functions of them. It may help to use categorical_logit_lpmf in place of categorical_lpmf.
  • Truncate one mixture component at $\mu$ from above and the other at $\mu$ from below, depending on your perspective on the dilemma raised by Jarle Tufto in the comments. It seems like you could get VERY different estimates depending on how exactly you decide to handle this. A nice sanity check: generate a fake dataset from the fitted parameters and make sure it has the right amount at 0, amount above $\mu$, etc.

Once you have a file with the right STAN code, you can use STAN with lots of different toolchains. To use it with R, check out these examples. I simplified one to get an MLE, using rstan::optimizing instead of sampling:

model = stan_model("Example1.stan")
fit = optimizing(model)

There are also some tricks for faster/better optimization that could help in practice.

  • 1
    $\begingroup$ If one uses optimizing, what is the underlying optimization that is performed? Is it something similar to gradient descent as in 'optim'? $\endgroup$ Commented Dec 27, 2021 at 19:14
  • 1
    $\begingroup$ It's L-BFGS by default. BFGS and Newton's method are available. $\endgroup$ Commented Dec 28, 2021 at 2:43

My STANswer is so complex that it's just begging for something to go wrong. Here's a simpler way: do all of your inference conditional on the (known) facts of whether each datum exceeds 0 and whether each datum exceeds $\mu$. In other words, reduce the data to:

  • The set of observations $S_1 \equiv \{y_i: 0<x_i<\mu\}$.
  • The set of observations $S_2 \equiv \{y_i: \mu<x_i\}$.
  • The number of zeroes $N_0$.
  • Let $N_1 \equiv |S_1|$ and $N_2 \equiv |S_2|$.

Then maximize, separately:

  • $(1-F_A(\mu))^{N_1} \prod_{y \in S_1} f_A(y) $ w.r.t. parameters of $f_A$.
  • $\prod_{y \in S_2} GPD(y)$ w.r.t. parameters of your gross promestic doduct (GPD).
  • $\pi^{N_0}(1-\pi)^{N_1 + N_2}$ w.r.t. $\pi$ .

This doesn't seem to do quite what you ask because:

  • $\mu$ must be user-specified.
  • the GPD scale parameter is not fixed to $\frac{f_A(\mu)}{1-F_A(\mu)}$. Hopefully that part is not essential for interpretability. If it is, maybe it would be good enough to just fix $\sigma$ based on the results of bullet point 1, then optimize the remaining parameters. It's no longer a joint MLE then. There's no way to uncouple the optimization if you are really dead set on that.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.