I have a Multivariate Normal Distribution with the mean vector $\boldsymbol\mu\in \mathbb{R}^n$ and the covariance matrix $\boldsymbol\Sigma \in \mathbb{R}^{n\times n}$ given as $$ f(\boldsymbol{x},\boldsymbol\mu,\boldsymbol\Sigma) = \det(2\pi\boldsymbol\Sigma)^{-\frac{1}{2}} \, \exp \left( -\frac{1}{2}(\boldsymbol{x} - \boldsymbol\mu)^{{{\!\mathsf{T}}}} \boldsymbol\Sigma^{-1}(\boldsymbol{x} - \boldsymbol\mu) \right) $$ Now, I want to compute the probability that a realization of $\boldsymbol{x}$ lies in a given polytopic set of the form $$ \mathcal{P} = \left\{ \boldsymbol{x} \;\middle|\; \boldsymbol{A} \boldsymbol{x} < \boldsymbol{b} \right\} $$ where the matrix $\boldsymbol{A} \in \mathbb{R}^{n\times m} $ and the vector $\boldsymbol{b} \in \mathbb{R}^{m} $ describes $m$ half-spaces and therefore a convex set.
The probability which I want to compute is given as $$ \Pr(\boldsymbol{A} \boldsymbol{x} < \boldsymbol{b}) = \int_\mathcal{P} f(\boldsymbol{x},\boldsymbol\mu,\boldsymbol\Sigma) \mathrm{d}\boldsymbol{x} $$
How can I compute/approximate this integration numerically in MATLAB or Python for a given mean $\boldsymbol\mu$, a given covariance $\boldsymbol\Sigma$, and a given set $\mathcal{P}$. In my problem, the variable $\boldsymbol{x}$ has around 10 to 100 dimensions.