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}(y;\mu,\sigma,\xi,\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?