Compute Probability of a Multivariate Normal Distribution over Polytope

I have a Multivariate Normal Distribution with the mean vector and the covariance matrix given as

Now, I want to compute the probability that a realization of lies in a given polytopic set of the form

where the matrix and the vector describes m half-spaces and therefore a convex set.

The probability which I want to compute is given as

How can I compute/approximate this integration numerically in MATLAB for a given mean μ, a given covariance Σ, and a given set 𝒫. In my problem, the variable x has around 10 to 100 dimensions.

Edit: ->

[Torsten]

Can the polytope be unbounded ? Or is it contained in a bounded box ?

[Michael Fink]

The polytope is bounded in my case.

[Torsten]

It could be quite time-consuming, but maybe Monte-Carlo integration is an option:

https://en.wikipedia.org/wiki/Monte_Carlo_integration

[Michael Fink]

Thank you @Torsten for the idea

Answer 1

[Paul]

Hi Michael,

If A is nonsingular, perhaps a change of coordinates will work

% z = A*x
muz = A*mux;
Sigmaz = A*Sigmax*A.';
ProbAxLTb = mvncdf(-inf+b,b,muz,Sigmaz);

See the doc page for mvncdf for examples, info, options, etc.

[Michael Fink]

Thanks Paul for the suggestion. However, my matrix A is generally not square. Typically it is a tall matrix. Therefore Sigmaz would be singular.

[Paul]

My mistake. I didn't squint hard enough to see A isn't square. But shouldn't A be m x n, not n x m as in the Question?

[Michael Fink]

Thanks Paul, yes there was also a mistake from my side. I fixed it in the first post.

[Paul]

Just want to make sure I'm clear on the question. For s simple case we have:

rng(100);
Sigma = rand(2);Sigma = Sigma*Sigma.';
mu = [1 1];
A = rand(4,2); b = (1:4).';
N = 1e5;
x = mvnrnd(mu,Sigma,N);
z = A*x.';
P = sum(all(z < b,1))/N

P = 0.8825

Is that the probabilty that you're trying to get through integration or other means?

[Michael Fink]

Hello Paul,

yes, exactly, this is the probability. But in my case, the Polytope (i.e. matrix A and vector b) is given/is deterministic and not random.

For example an unitbox for n=2 (Vertecise: (1,1),(-1,1),(1,-1),(-1,1) ) would be given as

A = [eye(2);-eye(2)]

A = 4×2

1 0 0 1 -1 0 0 -1

b = ones(4,1)

b = 4×1

1 1 1 1

[Torsten]

It was only an example.

Of course, you can directly use your values for sigma, mu, A and b in Paul's code instead of the randomly created.