# 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?

• May 2, 2016 at 11:22
• Oh. No, I had not. May 3, 2016 at 4:22
• On my computer, rHuber seems to work. What do you mean by 'rubbish'? Does it return letters? May 4, 2016 at 12:16
• @user603 I've done some tests and got unsatisfactory results but now I see it seems to work...
– Tim
May 7, 2016 at 20:33

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