1
$\begingroup$

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!

$\endgroup$
5
  • $\begingroup$ I'm seeing the same issue in V12.2 $\endgroup$ – mikado May 4 at 19:08
  • $\begingroup$ @mikado, I'm glad to hear its not just me . . . except that this then means Mathematica is broken(?) $\endgroup$ – David May 4 at 19:53
  • 2
    $\begingroup$ 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]. $\endgroup$ – JimB May 4 at 19:59
  • $\begingroup$ @JimB, yes this works! $\endgroup$ – David May 4 at 20:10
  • $\begingroup$ 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]. $\endgroup$ – JimB May 4 at 20:18
0
$\begingroup$

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

$\endgroup$
0
$\begingroup$

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}
]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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