6
$\begingroup$

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?

$\endgroup$
4
  • 1
    $\begingroup$ inside-r.org/packages/cran/marg/docs/pHuber ? $\endgroup$ May 2, 2016 at 11:22
  • $\begingroup$ Oh. No, I had not. $\endgroup$ May 3, 2016 at 4:22
  • $\begingroup$ On my computer, rHuber seems to work. What do you mean by 'rubbish'? Does it return letters? $\endgroup$
    – user603
    May 4, 2016 at 12:16
  • $\begingroup$ @user603 I've done some tests and got unsatisfactory results but now I see it seems to work... $\endgroup$
    – Tim
    May 7, 2016 at 20:33

1 Answer 1

6
$\begingroup$

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: enter image description here

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)))}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.