Simulate skew normal with given mean, variance & skewness

Question

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?

Answer 1

Implemented what Xi'an suggested. Well, at least i feel i learned something about this pdf.

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

}