Generating random samples from Huber density Huber density, connected to Huber loss, can be defined as:
$$ f(x) = \frac{1-\epsilon}{\sqrt{2\pi}} e^{-\rho_k(x)} $$
where
$$ \rho_k(x) =
   \begin{cases}
   \frac{1}{2} x^2       & |x|\le k \\
   k|x|- \frac{1}{2} k^2 & |x|>k
   \end{cases} $$
and $\epsilon$ satisfies
$$ \frac{2\phi(k)}{k} - 2\Phi(-k) = \frac{\epsilon}{1-\epsilon} $$
What is the best way to generate random samples for such distribution? Can we anyhow use the fact that it is kind of normal-Laplace mixture?
 A: As you suggest, this distribution is a mixture of a truncated Normal distribution and of a truncated Laplace distribution: namely,
$$f(x)\propto \exp\left\{-x^2/2\right\}\mathbb{I}_{(-k,k)}(x)+\exp\left\{-k|x|+k^2/2\right\}\mathbb{I}_{(-k,k)^c}(x)$$ implies that the distribution is the mixture of the Normal distribution truncated to $(-k,k)$ and of the Laplace distribution with rate $k$ truncated to $(-k,k)^c$. The only missing item is the weight $\alpha$ of the truncated Normal, which amounts to normalise both terms:
$$\exp\left\{-x^2/2\right\}=\frac{\sqrt{2\pi}}{\sqrt{2\pi}}\times
\dfrac{\Phi(k)-\Phi(-k)}{\Phi(k)-\Phi(-k)}\times\exp\left\{-x^2/2\right\}$$ hence the coefficient of the truncated Gaussian is
$$\sqrt{2\pi}\times\{\Phi(k)-\Phi(-k)\}$$while
$$\exp\left\{-k|x|+k^2/2\right\}=e^{k^2/2}\times\dfrac{2k^{-1}e^{-k^2}}{2k^{-1}e^{-k^2}}\times\exp\left\{-k|x|\right\}$$
since
$$\int_k^\infty \exp\left\{-k|x|\right\}\text{d}x=k^{-1}e^{-k^2}$$
Therefore this distribution is equal to
$$\alpha\mathcal{N}_{(-k,k)}(0,1)+(1-\alpha)\mathcal{L}_{(-k,k)^c}(0,k)$$with
$$\alpha=\dfrac{\sqrt{2\pi}\{\Phi(k)-\Phi(-k)\}}{\sqrt{2\pi}\{\Phi(k)-\Phi(-k)\}+2k^{-1}e^{-k^2/2}}$$
Simulating one of the truncated distributions is straightforward by cdf inversion. Here is an illustration of the fit for k=1:

based on the R code
genhuber <- function(k=1,n=1){
pk=pnorm(k);pmk=pnorm(-k);dk=pk-pmk
alp=1/(1+2*dnorm(k)/(k*dk))
u=runif(n)
return((u<alp)*qnorm(pmk+runif(n)*dk)+(u>alp)*(1-2*(runif(n)<.5))*(k+rexp(n)/k))
} 

dhuber <-function(x,k=1){
x=abs(x)
meps=1/(1+2*dnorm(k)/k-2*pnorm(-k))
return(meps*((x<k)*dnorm(x)+(x>=k)*exp(-k*(x-.5*k))/sqrt(2*pi)))}

