There is BinormalDistribution that handles easily bivariate normal distribution and its mixtures. There is also SkewNormalDistribution (univariate) that works nicely. But there's no such thing like BivariateSkewNormalDistribution, and I want to define it.

Its mathematical definition is


where $\varphi_p$ is the PDF of the $p$-variate normal distribution with mean $\mu$ and covariance matrix $\Sigma$; $\Phi_1$ denotes the CDF of a standard Gaussian; $y,\mu,a$ are 2-dimensional vectors. This looks like a formula that is straightforward to implement, so I did

Σ = {{σ1^2, ρ σ1 σ2}, {ρ σ1 σ2, σ2^2}};
s = Simplify@Sqrt[Det[Inverse[Σ]]];
t = Simplify@Sqrt[Tr[Inverse[Σ]] + 2 s];
sqrtΣ = Simplify[1/t {{Σ[[1, 1]] + s, Σ[[1, 2]]}, {Σ[[2, 1]], Σ[[2, 2]] + s}}];

SN[μ1_, μ2_, σ1_, σ2_, a1_, a2_, ρ_] := 
    PDF[BinormalDistribution[{μ1, μ2}, {σ1, σ2}, ρ], {x, y}]*
    (1 + Erf[{a1,a2}.sqrtΣ.({x, y} - {μ1, μ2})]),
    {x, -Infinity, Infinity}, {y, -Infinity, Infinity},
    Assumptions -> {σ1 > 0, σ2 > 0}]

I don't see a reason why this shouldn't work, but neither of the following works:

RandomVariate[ SN[μ1a, μ2a, σ1a, σ2a, a1a, a2a, ρa] /. {μ1a -> -0.233, μ2a -> 0.74, σ1a ->
 0.541, σ2a -> 0.259, a1a -> 0.7, a2a -> 0.5, ρa -> 0.049}]

Plot3D[ Evaluate@  PDF[SN[μ1a, μ2a, σ1a, σ2a, a1a, 
 a2a, ρa] /. {μ1a -> -0.233, μ2a -> 0.74, σ1a -> 0.541, σ2a -> 0.259, a1a -> 0.7, 
 a2a -> 0.5, ρa -> 0.049}, {x, y}], {x, -2., 2.}, {y, -1., 1.}]

I also put the full formula for sqrtΣ into the definition of SN, but I got an error message:

RandomVariate::noimp: Sampling from ProbabilityDistribution[1.13722 E^(-1.00241 (1.70834 Plus[<<2>>]^2-0.349703 (0.233 +\[FormalX]1) (-0.74+\[FormalX]2)+7.45368 Plus[<<2>>]^2)) (1+Erf[0.911681 (0.233 +\[FormalX]1)+0.631874 (-0.74+\[FormalX]2)]),{\[FormalX]1,-\[Infinity],\[Infinity]},{\[FormalX]2,-\[Infinity],\[Infinity]}] is not implemented. >>

A related question: Mixture of a bivariate distribution defined as a copula

  • $\begingroup$ Sampling from user defined multivariate distributions isn't implemented. You will need to somehow use built in distributions or define your own method for sampling. $\endgroup$
    – Andy Ross
    Jul 30, 2015 at 18:04
  • $\begingroup$ In general it is. I defined a sinh-arcsinh distribution (univariate though) and it works with no problems. Also, a mixture of binormal distributions works (even with the FingDistributionParameters, which I also want to perform for a mixture of binormal skew distributions as defined in my post). $\endgroup$
    – corey979
    Jul 30, 2015 at 18:09
  • 1
    $\begingroup$ Univariate is generally easy because one can use the inverse cdf if all else fails. Multivariate will probably never be implemented in full generality. Also, pretty much everything exists in closed form for multivariate normals so special case code exists for those. $\endgroup$
    – Andy Ross
    Jul 30, 2015 at 18:14
  • 1
    $\begingroup$ The Plot3D example doesn't work because the probability distribution needs to be pre-evaluated. Try adding an Evaluate in your code for SN. As for the random variates, you are still probably going to have to create your own method. $\endgroup$
    – Andy Ross
    Jul 30, 2015 at 18:32
  • 2
    $\begingroup$ It looks like your definition is a special case of "A multivariate skew normal distribution" (sciencedirect.com/science/article/pii/S0047259X03001313). And biomet.oxfordjournals.org/content/83/4/715.full.pdf+html pre-dates that article. $\endgroup$
    – JimB
    Nov 12, 2016 at 22:05

1 Answer 1


Maybe you might be able to find a bivariate density that matches what you want using CopulaDistribution:

cd = CopulaDistribution[{"Frank", 0.5}, {SkewNormalDistribution[0, 1, 0.2],
  SkewNormalDistribution[0, 1, 0.2]}];
PDF[cd, {x, y}]

$$\frac{0.0551589 \text{erfc}(-0.141421 x) \text{erfc}(-0.141421 y) e^{-\frac{x^2}{2}-\frac{y^2}{2}} 0.5^{-2 T(x,0.2)-2 T(y,0.2)+\frac{1}{2} \text{erfc}\left(-\frac{x}{\sqrt{2}}\right)+\frac{1}{2} \text{erfc}\left(-\frac{y}{\sqrt{2}}\right)}}{\left(1. 0.5^{-2 T(x,0.2)-2 T(y,0.2)+\frac{1}{2} \text{erfc}\left(-\frac{x}{\sqrt{2}}\right)+\frac{1}{2} \text{erfc}\left(-\frac{y}{\sqrt{2}}\right)}-1. 0.5^{\frac{1}{2} \text{erfc}\left(-\frac{x}{\sqrt{2}}\right)-2 T(x,0.2)}-1. 0.5^{\frac{1}{2} \text{erfc}\left(-\frac{y}{\sqrt{2}}\right)-2 T(y,0.2)}+0.5\, \right)^2}$$

RandomVariate[cd, 10]
(* {{3.27996, 0.610545}, {-0.162205, 0.789233}, {1.13898, -0.0868688},
   {0.283741, 0.396361}, {0.595174, -0.72346}, {0.110216, 1.39394}, 
   {1.52232,  0.0436556}, {-1.12374, 1.04785}, {-0.903327, -0.0190759},
   {0.730456, 0.66757}} *)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.