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

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

 - Large mass at zero
 - 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$:

$$\text{logit}(\pi)=f(\beta Z)$$

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 \\
\text{GPD}(y;\mu,\sigma,\xi),\,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.

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)$$

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

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

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