How to define a bivariate skew normal distribution?

Question

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

$$f(y)=2\varphi_p(y|\mu,\Sigma)\Phi_1(a^T\Sigma^{-1/2}(y-\mu))$$

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_, ρ_] := 
  ProbabilityDistribution[
    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

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}} *)