6
$\begingroup$

I want to draw samples from a skew normal distribution as part of a matlab project of mine. I already implemented the CDF and PDF of the distribution, but sampling from it still bothers me.
Sadly, the description of this process from the documentation of an R package is riddled with dead links, so I did some reading on the process.

One way of sampling from the distribution would be inverse transform sampling, which uses a uniform random variable $U\sim Unif(0,1) $ and involves solving

$F(F^{-1}(u)) = u$

with $F(x)$ being the CDF of the distribution we want to sample from. Since I don't know how to find the inverse of $F(x)$ myself, I did some searching, finding this question asked several times but not answered.

Edit: another method to sample from the distribution would be rejection sampling, however for that I need to find a distribution that

  1. I can draw samples from
  2. Has a pdf "which is at least as high at every point as the distribution we want to sample from, so that the former completely encloses the latter." (from the rejection sampling wiki article)

I've plotted the skew normal distribution with $\xi=1,\omega=1.5,\alpha=4$ and its truncated version (truncated to [0,2.5] here).

unrestricted and truncated skew-normal distribution

In my application of this, I will always truncate the distribution to a certain interval, so I'd need to find a distribution that 'contains' the SN pdf for (hopefully) all parameters.

Any ideas how I could go about sampling from such truncated skew-normal distributions?

$\endgroup$

migrated from math.stackexchange.com Nov 29 '17 at 16:41

This question came from our site for people studying math at any level and professionals in related fields.

  • $\begingroup$ What update are you referring to? I merely edited the tags to make the question better visible for others looking to sample from the SN distribution, that's why I replaced the normal with the skew normal tag. Your answer was what helped me resolve this, that's why I accepted it back in november. What's the problem here? $\endgroup$ – Rickyfox Jul 27 '18 at 12:46
  • $\begingroup$ OK, apologies then, I though you had moved from skew Normal to truncated skew Normal... $\endgroup$ – Xi'an Jul 27 '18 at 16:33
8
$\begingroup$

The most direct way of simulating a random variable from a distribution with cdf $F$ is to first simulate a Uniform variate $U\sim\mathcal{U}(0,1)$ and second return the inverse cdf transform $F^{-1}(U)$. When the inverse $F^{-1}$ is not available in closed form, a numerical inversion can be used. Numerical inversion may however be costly, especially in the tails.

One can also use accept-reject algorithms when the density $f$ is available and dominated by another density $g$, i.e., that there exists a constant $M$ such that$$f(x)<M g(x)$$. In the case of the skew-Normal distribution, the density is$$f(x)=2\varphi(x)\Phi(\alpha x)$$when $\varphi$ and $\Phi$ are the pdf and cdf of the standard Normal distribution, respectively. (Adding a location and a scale parameter does not modify the algorithm, since the outcome simply needs to be rescaled and translated.)

This density seems ideally suited for accept-reject since$$2\varphi(x)\Phi(\alpha x)< 2\varphi(x)$$as $\Phi$ is a cdf. This inequality implies that a first option to run accept-reject is to pick the Normal pdf for $g$ and $M=2$. This works out, as shown by the following picture when $\alpha=-3$:

enter image description here

and leads to an algorithm of the kind

T=1e3 #number of simulations
x=NULL
while (length(x)<T){
   y=rnorm(2*T)
   x=c(x,y[runif(2*T)<pnorm(alpha*y)])}
x=x[1:T]

which returns a reasonable fit of the pdf by the histogram:

enter image description here

There are however transforms of standard distributions that result in skew Normal variates: If $X_1,X_2$ are iid $\mathcal{N}(0,1)$, then

  1. $$\dfrac{\alpha|X_1|+X_2}{\sqrt{1+\alpha^2}}$$
  2. $$\dfrac{1+\alpha}{\sqrt{2(1+\alpha^2)}}\max\{X_1,X_2\}+\dfrac{1-\alpha}{\sqrt{2(1+\alpha^2)}}\min\{X_1,X_2\}$$

are skew Normal variates with parameter $\alpha$. (Both representations are identical when considering that $(X_1+X_2,X_1-X_2)/\sqrt{2}$ is an iid $\mathcal{N}(0,1)$ pair.)

$\endgroup$
2
$\begingroup$

I stumbled over this via googling. Rejection sampling is not needed. Instead, it is sufficient to flip the sign if the sample would be rejected!

This is because we can use that $\Phi(-ax)+\Phi(ax)=1$ and thus

$(f(x)+f(-x))/2= \phi(x) \Phi(ax)+ \phi(-x) \Phi(-ax)= \phi(x)$

Therefore, we can sample a skew-normal random variable by first sampling a standard normal variable $x$ and than flip its sign with probability $1-\Phi(ax)=\Phi(-ax)$. Beautiful!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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