Simulating from a skew normal distribution

Question

I am analyzing a data set in R, the histogram gives an impression of a normal distribution, but the qqplot suggests a slightly skewed normal distribution, so I want to try this out as well. I fitted the model like so:

library(sn)# Azzalini skewnormal package
m1  <-  selm(log(dataset) ~ 1,data=Data)

Is this the right way to do it? I also want to simulate data from this model for control, how do I do this? The "simulate" command does not work here.

Answer 1

The code you show is from the Azzalini skew normal (sn) package. To simulate data from the estimated model you can get estimated predicted values from the model, and then simulate from the distribution of the residuals, and add that. Since you didn't post your data I simulate some:

set.seed(7*11*13) 
test  <-  rsn(1000, 0, 1, alpha=1)
hist(test)
mod  <-  selm(test ~ 1)
mod
Object class: selm 
Call: selm(formula = test ~ 1)
Number of observations: 1000 
Number of covariates: 1 (includes constant term)
Number of parameters: 3 
Family: SN 
Estimation method: MLE
Log-likelihood: -1196.029 
summary(mod)
Call: selm(formula = test ~ 1)
Number of observations: 1000 
Family: SN 
Estimation method: MLE
Log-likelihood: -1196.029 
Parameter type: CP 

CP residuals:
     Min       1Q   Median       3Q      Max 
-2.66726 -0.56187 -0.07027  0.51962  2.68867 

Regression coefficients
     estimate  std.err  z-ratio Pr{>|z|}
mean  0.58863  0.02545 23.13155        0

Parameters of the SEC random component
       estimate std.err
s.d.     0.8050   0.019
gamma1   0.2746   0.080

Now you can get fitted values from this model by:

mod.fit  <-  fitted(mod)
str(mod.fit)
 Named num [1:1000] 0.589 0.589 0.589 0.589 0.589 ...
 - attr(*, "names")= chr [1:1000] "1" "2" "3" "4" ...

you can see all of them is equal to the fitted mean. Now you could continue, simulate new residuals from the fitted skew normal error distribution, and simulate from that. For the details you would have to look up Azzalini' book.

Answer 2

There is a distribution known as the skew normal distribution, with density $$f(x|\alpha,\omega,\xi)=\frac{2}{\omega}\,\varphi(\{x-\xi\}/\omega)\,\Phi(\alpha\{x-\xi\}/\omega)\qquad > x,\xi,\alpha\in\mathbb{R}\quad\omega\in\mathbb{R}^*_+$$ which generalises the Normal $\mathcal{N}(\xi,\omega)$ distribution when $\alpha\ne 0$.

You may try to fit this distribution to your dataset and test whether or not $\alpha=0$.