title: multivariate skew normal in R
body: <p>I'm trying to generate random numbers with a multivariate skew normal distribution using the <code>rmsn</code> command from the <code>sn</code> package in R. I would like, ideally, to be able to get three columns of numbers with a specified variances and covariances, while having one column strongly skewed. But I'm struggling to achieve both goals simultaneously.</p>
  
  <p>The post at <a href="https://stackoverflow.com/questions/4656784/skew-normal-distribution">skew normal distribution</a> was related and useful (and the source of some of the code below), but hasn't completely clarified the issue for me.</p>
  
  <p>I've been trying:</p>
  
  <pre><code>a &lt;- c(5, 0, 0) # set shape parameter
  s &lt;- diag(3) # create variance-covariance matrix
  w &lt;- sqrt(1/(1-((2*(a^2)/(1 + a^2))/pi))) # determine scale parameter to get sd of 1
  xi &lt;- w*a/sqrt(1 + a^2)*sqrt(2/pi) # determine location parameter to get mean of 0
  
  apply(rmsn(n=1000, xi=c(xi), Omega=s, alpha=a), 2, sd)
  colMeans(rmsn(n=1000, xi=c(xi), Omega=s, alpha=a))
  </code></pre>
  
  <p>The columns means and SDs are correct for the second and third columns (which have no skew) but not the first (which does). Can anyone clarify where my code above, or my thinking, has gone wrong? I may be misunderstanding how to use <code>rmsn</code>, or the output. Any assistance would be appreciated.</p>