How to fit a skew normal distributon to given data? I've got some data which I want to fit it to a skew normal distribution given by
$f(z)=\frac{2}{\sigma}\phi(\frac{z-\mu}{\sigma})\Phi(\lambda\frac{z-\mu}{\sigma})$
where $\phi(z)=\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}z^2}$,$\Phi(z)=\int_{-\infty}^{z}\phi(t)dt$ and $\lambda$ is a factor that control skewness.
Question is: how should i do to fit my data to such a formula?Other words ,which method can i use to estimate the $\mu,\sigma$ and $\lambda$?
Any insights or suggestions will be appreciated!
 A: You can use MLE (maximum likelihood estimation), see also the links in comments, where there are given some formulas for method of moments. The easy way out is to use the sn package in R (on CRAN).  sn use MLE (maximum likelihood estimation).
Since you didn't give any data or context, I will just simulate some data, and show the use of the package:
    library(sn)
    set.seed(7*11*13) #My public seed
    testdata <- data.frame(X=rsn(100, omega=3, alpha=0.5))
    mod <- selm(X ~ 1, data=testdata)
    
    summary(mod)
    Call: selm(formula = X ~ 1, data = testdata)
    Number of observations: 100 
    Family: SN 
    Estimation method: MLE
    Log-likelihood: -249.2581 
    Parameter type: CP 
    
    CP residuals:
         Min       1Q   Median       3Q      Max 
    -6.54250 -2.11290 -0.01777  2.18490  8.39505 
    
    Regression coefficients
         estimate std.err z-ratio Pr{>|z|}
    mean   1.7614  0.2927  6.0170        0
    
    Parameters of the SEC random component
           estimate std.err
    s.d.    2.92717   0.208
    gamma1  0.06278   0.257  

