Multivariate normal distribution using python scipy stats and integrate nquad

Let

be independent normal random variables with means

and unit variances, i.e.

I would like to compute the probability

Is there any easy way to compute the above probability using scipy.stats.multivariate_normal? If not, how do we do it using scipy.integrate?

Based on the info from one of SO post mentioned in the comments: https://math.stackexchange.com/questions/270745/compute-probability-of-a-particular-ordering-of-normal-random-variables

Also same method mentioned on wiki here: https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Affine_transformation

You should convert F(X) to F(Y) in the following way:

import numpy as np
from scipy.stats import norm, multivariate_normal


mvn = multivariate_normal



def compute_probability(thetas):
    """
        following:
        https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Affine_transformation


        F(X1<X2<X3<...<Xn) = F(Y1<0, Y2<0, ... Y_{n-1}<0)

        where:
            Y_i = X_{i-1} - X_{i}
            ...
        
        or what is the same:
            Y = c + BX
            mean and sigma of Y can be found based on formulas from wikipedia

    """
    
    n = len(thetas)

    # set diagonal to 1
    B = np.eye(n)

    # set right to diagonal to -1
    idx = np.arange(n-1)
    B[idx,idx+1] = -1

    # remove last row
    B = B[:-1]

    # calculate multivate sigma and mu
    sigma = np.eye(n)
    mu_new = B.dot(thetas.T)


    sigma_new = B.dot(sigma).dot(B.T)


    MVN = mvn(mean=mu_new, cov = sigma_new)

    x0 = np.zeros(shape = n-1)
    p = MVN.cdf(x = x0)
    return p


# Example usage:
theta = np.array([1, -2, 0.5])  # Example coefficient
p = compute_probability(theta)

print(p)

Outputs:

theta = (0,0)
p = 0.5

theta = (0,0,0)
p = 0.1666

theta = (100, -100)
p = 0

theta = (-100, 100)
p = 1

theta = (0,0,0,100, -100)
p = 0