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

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!

• I'm seeing the same issue in V12.2 – mikado May 4 at 19:08
• @mikado, I'm glad to hear its not just me . . . except that this then means Mathematica is broken(?) – David May 4 at 19:53
• Not broken. There is no nice and compact answer although Mathematica keeps trying. The tmvtnorm library does it numerically which is a big difference. Try NExpectation[{x, y}, {x, y} \[Distributed] dTruncated]. – JimB May 4 at 19:59
• @JimB, yes this works! – David May 4 at 20:10
• Mathematica doesn't use the same parameter structure so your examples don't match up on the bivariate normal or the truncation limits. To match the R code you'd need the following: d = BinormalDistribution[{0, 0}, {0.5^0.5, 1}, 0.7071067811865475] or d = MultinormalDistribution[{0, 0}, {{0.5, 0.5}, {0.5, 1}}] and dTruncated = TruncatedDistribution[{{-0.5, 0}, {-\[Infinity], 2}}, d]. – JimB May 4 at 20:18

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

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}
]
`