h(T_j) = x_j^{\top}\bm{\beta} + \varepsilon_j, \,\,\, j=1,\dots,n,\end{aligned}$$

where ${\bf X}= (x_1,\dots,x_n)^{\top}$ is an $n\times p$ known design matrix and $\varepsilon_j \stackrel{ind.}{\sim} F(\cdot;\bm{\theta})$, $F$ is a continuous distribution with support on ${\mathbb R}$ and parameters $\bm{\theta}\in \Theta \subset {\mathbb R}^d$, and $h:{\mathbb R}_+\rightarrow{\mathbb R}$ is a continuous increasing function. The most common choice for $h$ is the logarithmic function, while the distribution of the errors $\varepsilon_j$ is typically assumed to be normal. Given that the assumption of normality of the errors can be restrictive in practice, other distributional assumptions have been recently studied such as the log Birnbaum-Saunders [@B08], finite mixtures of normal distributions [@KL08], the symmetric family of scale mixtures of normals [@VS15], and the log-generalised extreme value distribution [@RD15].

AFT models are extremely relevant in medicine, given that survival data naturally arise in many medical studies, which typically involve the follow-up of other covariates. The presence of different types of censored observations is common in this context [@B08; @KL08; @VS15]. If we assume that the errors $\varepsilon_j$ are distributed according to a LTP distribution with location $0$ and $\bm{\theta}=(\sigma,\gamma,\delta)$, then we can write the likelihood function as follows, $$\begin{aligned}
L(\bm{\beta},\sigma,\gamma,\delta) =& \prod_{j\in\text{Obs}} s_l(T_j;x_j^{\top}\bm{\beta},\sigma,\gamma,\delta) \times \prod_{j\in\text{Left}} S_l(T_j;x_j^{\top}\bm{\beta},\sigma,\gamma,\delta)\\ &\times  \prod_{j\in\text{Right}}\left[1-S_l(T_j;x_j^{\top}\bm{\beta},\sigma,\gamma,\delta)\right]\\
&\times  