# Simulate skew normal with given mean, variance & skewness

In short, I would like to find the parameters of the skew normal so that given first three moments are matched. Find location, scale and shape of skew normal, so that mean = 0, variance = 1 and skewness = arbitrary value.

I would like to simulate from a skew-normal with the following properties

$E[Y] = 0 \\ E[Y^2] = 1\\ E[Y^3] = a \qquad a \quad \text{being any (possible) real number}$

I would say this is the same as a r.v. with mean = 0, variance = 1 and skewness = a. I am trying to use a skew normal Source1, wikipedia

There are formulas for the mean, variance and skewness. My approach was solving this with nleqlsv in R, but my optimization stops far away from a good result.

This is my whole R Code

library(sn)

sn_mean <- function(loc, scale, shape) {
delta <- shape / sqrt(1 + shape^2)
return ( loc + scale * delta * sqrt(2/ pi) )
}

sn_var <- function(loc, scale, shape) {
delta <- shape / sqrt(1 + shape^2)
return ( scale^2 * (1 - 2 * delta^2 / pi) )
}

sn_skew <- function(loc, scale, shape){
delta <- shape / sqrt(1 + shape^2)
return ( (4-pi) / 2 ) * ( ( delta * sqrt(2/pi))^3 / (1 - 2* delta^2 / pi)^(3/2) )
}

skew <- function(X){
return (sum(1/length(X) * ((X-mean(X))/sqrt(var(X)))^3))
}

objfunc <- function(x) {
shape <- x[1]
scale <- x[2]
loc <- x[3]
y <- numeric(3)
y[1] <- 0 - sn_mean(loc, scale, shape)
y[2] <- 1 - sn_var(loc, scale, shape)
y[3] <- 1 - sn_skew(loc, scale, shape)
return(y)
}

a <- nleqslv(c(100,100,100), objfunc,  control=list(allowSingular=TRUE))

• Is there a better approach to get a skew normal with specific mean, var and skew?
• How can I improve my optimization to yield a better result?
• Possible duplicate of Sampling from Skew Normal Distribution – kjetil b halvorsen Dec 1 '17 at 18:41
• I don't think this is a duplicate. There is a already a function for simulating from skew normal in r. I am looking to get the right parameters for such a function not how to simulate if I have the parameters. Mine is about parameter fitting based on first three moments. – PalimPalim Dec 1 '17 at 18:43
• @kjetilbhalvorsen please check my edits and if you agree remove duplicate tag. – PalimPalim Dec 1 '17 at 18:47
• Since a skew Normal has three parameters and you have three constraints, where is the issue? – Xi'an Dec 1 '17 at 21:30
• Skewness only depends on $\delta$, hence identifies $\delta$ if there exists a solution. Then variance only depends on $\omega$ and $\delta$, hence identifies $\omega$ if there exists a solution. And mean therefore identifies $\mu$. (I use the notations of Wikipedia.) – Xi'an Dec 1 '17 at 22:02

There are a few implicit bounds for each parameter. PalimPalim tried to solve a <- nleqslv(c(100,100,100), objfunc, control=list(allowSingular=TRUE)), but skewness=100 is prohibited. Skewness is bounded to be -.9952717 to +.9952717, approximately, $$\frac{4-\pi}{2}\left (\frac{2}{\pi-2} \right )^{3/2}$$ to be exact. And negative variance does not make sense either.

The function below gives you back the three parameter for skew normal, given a vector of length 3, having mean, variance and skewness that you wish.

invsn <- function(mvs) {
# given a vector of mean, variance, skenewss, returns vector of xi, omeaga, alpha
# of skew normal distribution.
mu <- mvs[1]
var <- mvs[2]
sk <- mvs[3]

sn_skew <- function(delta_bounded) {
# delta has to stay in (-.5*pi, .5*pi) range.
# i came up with use of atan() to bound the value
delta <- sqrt(abs(atan(delta_bounded)))
if (delta_bounded<0)  delta <- -delta
sk <- (4 - pi )/2 * (delta * sqrt(2/pi))**3 / ((1-2*(delta**2)/pi)**1.5)
sk
}
sk0 <- max(min(sk, .5), -.5) # arbitrary choice to not allow starting from extreme value
o <- nleqslv(sk0, function(x) sn_skew(x) - sk)
if (o$$termcd != 1) {print(o); stop('delta trouble')} delta_bounded <- o$$x
delta <- sqrt(abs(atan(delta_bounded)))
if (delta_bounded<0)  delta <- -delta
# definition of delta bounds it value to (-1,1): -alpha/sqrt(1+alpha**2) approaches
# -1 and 1 as alpha goes -Inf and +Inf. so...
if (delta <= -1 || delta >= 1) stop('delta out of bound:  skewness cant go beyond +/-.9952717')

# move on
if (var < 0) stop('negative variance doesnt make sense')
omega <- sqrt(var / ( 1 - 2*(delta**2)/pi))

xi <- mu - omega*delta*sqrt(2/pi)

sn_delta <- function(alpha) alpha/sqrt(1+alpha**2)
alpha0 <- delta
o <- nleqslv(alpha0, function(x) sn_delta(x) - delta)
if (o$$termcd != 1) {print(o); stop('alpha trouble')} alpha <- o$$x

ret <- c(xi, omega, alpha)
names(ret) <- c('xi', 'omega', 'alpha')
ret

}

• I tried the function for an example and confirm that the moments of the simulated sample are as expected. Great work! – PalimPalim Dec 19 '18 at 9:24