The Skew Normal Distribution

Let \(X_i \sim \text{SkewNormal}(\mu,\sigma,\lambda)\) (location (\({\mathbb R}\)), scale (\({\mathbb R}_+\)), and shape (\({\mathbb R}\)) parameters). The pdf of the skew Normal distribution is \[ s(x;\mu,\sigma,\lambda) = \dfrac{2}{\sigma}\phi\left( \dfrac{x-\mu}{\sigma} \right)\Phi\left(\lambda \dfrac{x-\mu}{\sigma}\right), \] Then, if \(\lambda = 0\), \(s=\phi\). Thus, we say that the Normal distribution is ``nested’’ in the skew Normal distribution. If \(X_i \sim \text{SkewNormal}(\mu,\sigma,\lambda)\), we may be interested on testing: \[H_0: \lambda = 0 \,\,\,\,\,\, vs. \,\,\,\,\,\,\, H_1: \lambda \neq 0,\] in order to decide whether the data are Normally or skew-Normally distributed.

The Likelihood Ratio Test is based on \(T_n = - 2 \log\lambda({\bf x})\), where \[\lambda({\bf x})= \dfrac{\sup_{\mu,\sigma} \prod_{j=1}^{n} \dfrac{1}{\sigma}\phi\left(\dfrac{x_j-\mu}{\sigma}\right)}{\sup_{\mu,\sigma,\lambda} \prod_{j=1}^{n} \dfrac{2}{\sigma}\phi\left(\dfrac{x_j-\mu}{\sigma}\right) \Phi\left(\lambda\dfrac{x_j-\mu}{\sigma}\right)}\] We reject the null hypothesis if \(T_n \geq \chi^2_{1,1-\alpha}\). The calculation of \(T_n\) requires numerical methods. Next, we present several examples to illustrate the use of this test for testing normality vs skew-normality.

Data generated from a SkewNormal(0,1,5)

rm(list=ls())
# Required packages
library(sn)
## Loading required package: stats4
## 
## Attaching package: 'sn'
## The following object is masked from 'package:stats':
## 
##     sd
library(knitr)

# Simulated data
set.seed(1234)
data.sim <- rsn(500,0,1,5)

hist(data.sim, breaks = 20, col = "gray")
box()

# Log-likelihood for the skew normal 
llsn <- function(par){
  if(par[2]>0) return( sum(dsn(data.sim, par[1], par[2], par[3], log = T))  )
  else return(-Inf)
}

# MLE for the skew normal parameters
# Three initial points
init1 <- c(mean(data.sim), sd(data.sim),-0.5)
init2 <- c(mean(data.sim), sd(data.sim),0)
init3 <- c(mean(data.sim), sd(data.sim),0.5)

OPT1 <- optim(init1, llsn, control = list(fnscale = -1, maxit = 10000))
OPT2 <- optim(init2, llsn, control = list(fnscale = -1, maxit = 10000))
OPT3 <- optim(init3, llsn, control = list(fnscale = -1, maxit = 10000))

indmax <- which.max(c(OPT1$value,OPT2$value,OPT3$value))
indmax # Best initial point (it produces the largest log likelihood)
## [1] 2
if(indmax==1) MLESN <- OPT1$par; if(indmax==2) MLESN <- OPT2$par; if(indmax==3) MLESN <- OPT3$par;
kable(MLESN, digits = 3)
x
-0.001
1.035
5.140 
# Log-likelihood for the normal 
lln <- function(par){
  if(par[2]>0) return( sum(dnorm(data.sim, par[1], par[2], log = T))  )
  else return(-Inf)
}

# MLE for the normal parameters
MLEN <- optim(c(0,1), lln, control = list(fnscale = -1))$par
kable(MLEN, digits = 3)
x
0.802
0.653 
# LR test statistic
LRT <- -2*(lln(MLEN) - llsn(MLESN))
kable(LRT, digits = 3)
x
96.09 
# 0.95 Quantile of a chi-square with 1 d.f.
q <- qchisq(0.95, df = 1)
kable(q, digits = 3)
x
3.841 
# Outcome from the test (approximate size = 0.05)
ifelse(LRT > q, "Null Hypothesis Rejected", "Null Hypothesis Not Rejected")
## [1] "Null Hypothesis Rejected"

Data generated from a Normal(0,1)

rm(list=ls())
# Required packages
library(sn)
library(knitr)

# Simulated data
set.seed(1234)
data.sim <- rnorm(500,0,1)

hist(data.sim, breaks = 20, col = "gray")
box()

# Log-likelihood for the skew normal 

llsn <- function(par){
  if(par[2]>0) return( sum(dsn(data.sim, par[1], par[2], par[3], log = T))  )
  else return(-Inf)
}

# MLE for the skew normal parameters
# Three initial points
init1 <- c(mean(data.sim), sd(data.sim),-0.5)
init2 <- c(mean(data.sim), sd(data.sim),0)
init3 <- c(mean(data.sim), sd(data.sim),0.5)

OPT1 <- optim(init1, llsn, control = list(fnscale = -1, maxit = 10000))
OPT2 <- optim(init2, llsn, control = list(fnscale = -1, maxit = 10000))
OPT3 <- optim(init3, llsn, control = list(fnscale = -1, maxit = 10000))

indmax <- which.max(c(OPT1$value,OPT2$value,OPT3$value))
indmax # Best initial point (it produces the largest log likelihood)
## [1] 3
if(indmax==1) MLESN <- OPT1$par; if(indmax==2) MLESN <- OPT2$par; if(indmax==3) MLESN <- OPT3$par;
kable(MLESN, digits = 3)
x
-0.550
1.172
0.732 
# Log-likelihood for the normal 
lln <- function(par){
  if(par[2]>0) return( sum(dnorm(data.sim, par[1], par[2], log = T))  )
  else return(-Inf)
}

# MLE for the normal parameters
MLEN <- optim(c(0,1), lln, control = list(fnscale = -1))$par
kable(MLEN, digits = 3)
x
0.002
1.034 
# LR test statistic
LRT <- -2*(lln(MLEN) - llsn(MLESN))
kable(LRT, digits = 5)
x
0.31952 
# 0.95 Quantile of a chi-square with 1 d.f.
q <- qchisq(0.95, df = 1)
kable(q, digits = 3)
x
3.841 
# Outcome from the test (approximate size = 0.05)
ifelse(LRT > q, "Null Hypothesis Rejected", "Null Hypothesis Not Rejected")
## [1] "Null Hypothesis Not Rejected"

Monte Carlo Simulation to check the size of the LRT (when the data are generated from a Normal(0,1))

rm(list=ls())
# Required packages
library(sn)
library(knitr)

indicator = vector()

for(i in 1:10000){
  # Simulated data
  set.seed(i)
  data.sim <- rnorm(500,0,1)
  
  # Log-likelihood for the skew normal 
  llsn <- function(par){
    if(par[2]>0) return( sum(dsn(data.sim, par[1], par[2], par[3], log = T))  )
    else return(-Inf)
  }
  
  # MLE for the skew normal parameters
  # Three initial points
  init1 <- c(mean(data.sim), sd(data.sim),-0.5)
  init2 <- c(mean(data.sim), sd(data.sim),0)
  init3 <- c(mean(data.sim), sd(data.sim),0.5)
  
  OPT1 <- optim(init1, llsn, control = list(fnscale = -1, maxit = 10000))
  OPT2 <- optim(init2, llsn, control = list(fnscale = -1, maxit = 10000))
  OPT3 <- optim(init3, llsn, control = list(fnscale = -1, maxit = 10000))
  
  indmax <- which.max(c(OPT1$value,OPT2$value,OPT3$value))
  indmax # Best initial point (it produces the largest log likelihood)
  
  if(indmax==1) MLESN <- OPT1$par; if(indmax==2) MLESN <- OPT2$par; if(indmax==3) MLESN <- OPT3$par;
  
  # Log-likelihood for the normal 
  lln <- function(par){
    if(par[2]>0) return( sum(dnorm(data.sim, par[1], par[2], log = T))  )
    else return(-Inf)
  }
  
  # MLE for the normal parameters
  MLEN <- optim(c(0,1), lln, control = list(fnscale = -1))$par
  
  # LR test statistic
  LRT <- -2*(lln(MLEN) - llsn(MLESN))
  
  # 0.95 Quantile of a chi-square with 1 d.f.
  q <- qchisq(0.95, df = 1)
  
  # Outcome from the test
  indicator[i] <- ifelse(LRT > q, 1, 0)
}

# Approximate size
mean(indicator)
## [1] 0.046

Monte Carlo Simulation to check the size of the LRT (when the data are generated from a SkewNormal(0,1,0.5))

## Monte Carlo Simulation to check the size of the LRT (when the data are generated from a SkewNormal(0,1,0.5)
rm(list=ls())
# Required packages
library(sn)
library(knitr)

indicator = vector()

for(i in 1:10000){
  # Simulated data
  set.seed(i)
  data.sim <- rsn(500,0,1,0.5)
  
  # Log-likelihood for the skew normal 
  llsn <- function(par){
    if(par[2]>0) return( sum(dsn(data.sim, par[1], par[2], par[3], log = T))  )
    else return(-Inf)
  }
  
  # MLE for the skew normal parameters
  # Three initial points
  init1 <- c(mean(data.sim), sd(data.sim),-0.5)
  init2 <- c(mean(data.sim), sd(data.sim),0)
  init3 <- c(mean(data.sim), sd(data.sim),0.5)
  
  OPT1 <- optim(init1, llsn, control = list(fnscale = -1, maxit = 10000))
  OPT2 <- optim(init2, llsn, control = list(fnscale = -1, maxit = 10000))
  OPT3 <- optim(init3, llsn, control = list(fnscale = -1, maxit = 10000))
  
  indmax <- which.max(c(OPT1$value,OPT2$value,OPT3$value))
  indmax # Best initial point (it produces the largest log likelihood)
  
  if(indmax==1) MLESN <- OPT1$par; if(indmax==2) MLESN <- OPT2$par; if(indmax==3) MLESN <- OPT3$par;
  
  # Log-likelihood for the normal 
  lln <- function(par){
    if(par[2]>0) return( sum(dnorm(data.sim, par[1], par[2], log = T))  )
    else return(-Inf)
  }
  
  # MLE for the normal parameters
  MLEN <- optim(c(0,1), lln, control = list(fnscale = -1))$par
  
  # LR test statistic
  LRT <- -2*(lln(MLEN) - llsn(MLESN))
  
  # 0.95 Quantile of a chi-square with 1 d.f.
  q <- qchisq(0.95, df = 1)
  
  # Outcome from the test
  indicator[i] <- ifelse(LRT > q, 1, 0)
}

# Approximate size
mean(indicator)
## [1] 0.0543

Low rejection rate when the null hypothesis is false.

Monte Carlo Simulation to check the size of the LRT (when the data are generated from a SkewNormal(0,1,1))

## Monte Carlo Simulation to check the size of the LRT (when the data are generated from a SkewNormal(0,1,1))
rm(list=ls())
# Required packages
library(sn)
library(knitr)

indicator = vector()

for(i in 1:10000){
  # Simulated data
  set.seed(i)
  data.sim <- rsn(500,0,1,1)
  
  # Log-likelihood for the skew normal 
  llsn <- function(par){
    if(par[2]>0) return( sum(dsn(data.sim, par[1], par[2], par[3], log = T))  )
    else return(-Inf)
  }
  
  # MLE for the skew normal parameters
  # Three initial points
  init1 <- c(mean(data.sim), sd(data.sim),-0.5)
  init2 <- c(mean(data.sim), sd(data.sim),0)
  init3 <- c(mean(data.sim), sd(data.sim),0.5)
  
  OPT1 <- optim(init1, llsn, control = list(fnscale = -1, maxit = 10000))
  OPT2 <- optim(init2, llsn, control = list(fnscale = -1, maxit = 10000))
  OPT3 <- optim(init3, llsn, control = list(fnscale = -1, maxit = 10000))
  
  indmax <- which.max(c(OPT1$value,OPT2$value,OPT3$value))
  indmax # Best initial point (it produces the largest log likelihood)
  
  if(indmax==1) MLESN <- OPT1$par; if(indmax==2) MLESN <- OPT2$par; if(indmax==3) MLESN <- OPT3$par;
  
  # Log-likelihood for the normal 
  lln <- function(par){
    if(par[2]>0) return( sum(dnorm(data.sim, par[1], par[2], log = T))  )
    else return(-Inf)
  }
  
  # MLE for the normal parameters
  MLEN <- optim(c(0,1), lln, control = list(fnscale = -1))$par
  
  # LR test statistic
  LRT <- -2*(lln(MLEN) - llsn(MLESN))
  
  # 0.95 Quantile of a chi-square with 1 d.f.
  q <- qchisq(0.95, df = 1)
  
  # Outcome from the test
  indicator[i] <- ifelse(LRT > q, 1, 0)
}

# Approximate size
mean(indicator)
## [1] 0.2331

Low rejection rate when the null hypothesis is false.

Monte Carlo Simulation to check the size of the LRT (when the data are generated from a SkewNormal(0,1,2))

## Monte Carlo Simulation to check the size of the LRT (when the data are generated from a SkewNormal(0,1,2))
rm(list=ls())
# Required packages
library(sn)
library(knitr)

indicator = vector()

for(i in 1:10000){
  # Simulated data
  set.seed(i)
  data.sim <- rsn(500,0,1,2)
  
  # Log-likelihood for the skew normal 
  llsn <- function(par){
    if(par[2]>0) return( sum(dsn(data.sim, par[1], par[2], par[3], log = T))  )
    else return(-Inf)
  }
  
  # MLE for the skew normal parameters
  # Three initial points
  init1 <- c(mean(data.sim), sd(data.sim),-0.5)
  init2 <- c(mean(data.sim), sd(data.sim),0)
  init3 <- c(mean(data.sim), sd(data.sim),0.5)
  
  OPT1 <- optim(init1, llsn, control = list(fnscale = -1, maxit = 10000))
  OPT2 <- optim(init2, llsn, control = list(fnscale = -1, maxit = 10000))
  OPT3 <- optim(init3, llsn, control = list(fnscale = -1, maxit = 10000))
  
  indmax <- which.max(c(OPT1$value,OPT2$value,OPT3$value))
  indmax # Best initial point (it produces the largest log likelihood)
  
  if(indmax==1) MLESN <- OPT1$par; if(indmax==2) MLESN <- OPT2$par; if(indmax==3) MLESN <- OPT3$par;
  
  # Log-likelihood for the normal 
  lln <- function(par){
    if(par[2]>0) return( sum(dnorm(data.sim, par[1], par[2], log = T))  )
    else return(-Inf)
  }
  
  # MLE for the normal parameters
  MLEN <- optim(c(0,1), lln, control = list(fnscale = -1))$par
  
  # LR test statistic
  LRT <- -2*(lln(MLEN) - llsn(MLESN))
  
  # 0.95 Quantile of a chi-square with 1 d.f.
  q <- qchisq(0.95, df = 1)
  
  # Outcome from the test
  indicator[i] <- ifelse(LRT > q, 1, 0)
}

# Approximate size
mean(indicator)
## [1] 0.9831