Fitting custom distributions by MLE

Question

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?

Answer 1

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);
    else
      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:

install.packages("rstan")
library("rstan")
model = stan_model("Example1.stan")
fit = optimizing(model)

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

Answer 2

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.