# Simulating from a skew normal distribution

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.

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.

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$.

• I don't think this answers the question ... Apr 27 '17 at 17:02