# Fitting custom distributions by MLE

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?

• 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. – Jarle Tufto Aug 25 '17 at 13:53
• 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. – Lucas Roberts Aug 25 '17 at 17:32
• @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)? – Lucas Roberts Sep 1 '17 at 13:13
• @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. – Lucas Roberts Sep 2 '17 at 12:59
• @StatsPlease -> thinking about this a bit more -> you could get at what you want via a Tweedie model (en.wikipedia.org/wiki/Tweedie_distribution) and that will give you the 0-point mass as well as the continuous on positive values distribution, you could then do a tail mixture approach that you want based on a GPD mixture. You could then add a zero-inflation component on top of the Tweedie but that may not be necessary once you've used the Tweedie. – Lucas Roberts Sep 20 '17 at 23:21