2
$\begingroup$

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

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?
$\endgroup$
  • 1
    $\begingroup$ Possible duplicate of Sampling from Skew Normal Distribution $\endgroup$ – kjetil b halvorsen Dec 1 '17 at 18:41
  • $\begingroup$ 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. $\endgroup$ – PalimPalim Dec 1 '17 at 18:43
  • $\begingroup$ @kjetilbhalvorsen please check my edits and if you agree remove duplicate tag. $\endgroup$ – PalimPalim Dec 1 '17 at 18:47
  • 2
    $\begingroup$ Since a skew Normal has three parameters and you have three constraints, where is the issue? $\endgroup$ – Xi'an Dec 1 '17 at 21:30
  • 2
    $\begingroup$ 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.) $\endgroup$ – Xi'an Dec 1 '17 at 22:02
1
$\begingroup$

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

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

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.