Mathematica can't seem to handle Truncated BinormalDistribution when there is non-zero correlation coefficient

Question

I would like to use Mathematica to analyze (e.g., compute moments, plot, etc) a truncated bivariate normal distribution. For example:

d = BinormalDistribution[{0,0},{.5,1},.5];
dTruncated = TruncatedDistribution[{{-.5,Infinity},{0,2}},d]
Mean[dTruncated]

When I run this code, though, Mathematica begins evaluating and never stops (I ran it all night and nothing). I don't get any error messages. Same when I try to plot the PDF of dTruncated or sample points from the distribution.

I'm running Mathematica v 11.2 with Windows 10.0 on a 4.6GHz Intel i9 processor with 64Gb RAM, so I don't think it's a processing speed issue.

The problem only seems to occur when the correlation coefficient is non-zero. When I run the same code as above but just make the correlation coefficient in BinormalDistribution = 0, it works fine:

d = BinormalDistribution[{0,0},{.5,1},0];
dTruncated = TruncatedDistribution[{{-.5,Infinity},{0,2}},d]
Mean[dTruncated]

This immediately spits out an answer. I have tried numerous combinations of parameter values, and it only ever works when the correlation coefficient equals 0. Unfortunately, that's not very helpful for me.

There is an R package that does this easily (see here and here) in a few lines of code:

> library(tmvtnorm)
> mu <- c(0, 0)
> sigma <- matrix(c(.5, .5, .5, 1), 2, 2)
> a <- c(-0.5, -Inf)
> b <- c(0, 2)
> moments <- mtmvnorm(mean=mu, sigma=sigma,
> lower=a, upper=b)

Any assistance with this would be very much appreciated!

Answer 1

There are workarounds:

d = BinormalDistribution[{0,0},{.5,1},.5];
dTruncated = TruncatedDistribution[{{-.5,Infinity},{0,2}},d]
Mean[MarginalDistribution[dTruncated,1]]

0.241775

 Mean[MarginalDistribution[dTruncated,2]]

Fail i.e. the returned input

and

PDF[dTruncated, {x, y}] // InputForm

Piecewise[ {{0.8238509428557775* E^(0.6666666666666666* (-4.*x^2 + 2.*x*y - y^2)), Inequality[-0.5, Less, x, LessEqual, Infinity] && Inequality[0, Less, y, LessEqual, 2]}}, 0]

NIntegrate[x*Piecewise[{{0.8238509428557775*
 E^(0.6666666666666666*(-4.*x^2 + 2.*x*y - y^2)), 
 Inequality[-0.5, Less, x, LessEqual, Infinity] && 
 Inequality[0, Less, y, LessEqual, 2]}}, 0], {x, -Infinity,Infinity}, {y, -Infinity, Infinity}]

0.241775

NIntegrate[y*Piecewise[{{0.8238509428557775*
 E^(0.6666666666666666*(-4.*x^2 + 2.*x*y - y^2)), 
 Inequality[-0.5, Less, x, LessEqual, Infinity] && 
 Inequality[0, Less, y, LessEqual, 2]}}, 0], {x, -Infinity,Infinity}, {y, -Infinity, Infinity}]

0.739385

Answer 2

This answer is only intended to show that the solution can be obtained by numerical integration is a single step. (Otherwise, this answer is essentially the same as the previous one.)

d = BinormalDistribution[{0,0},{.5,1},.5];
dTruncated = TruncatedDistribution[{{-.5,Infinity},{0,2}},d]

NIntegrate[
  {x, y} * PDF[dTruncated, {x, y}], 
  {x, -Infinity, Infinity},
  {y, -Infinity, Infinity}
]