Многомерное нормальное распределение с использованием статистики Python scipy и интеграции nquad

На основе информации из одного из сообщений SO, упомянутого в комментариях: https://math.stackexchange.com/questions/270745/compute-probability-of-a-private-ordering-of-normal-random-variables

Также тот же метод, упомянутый здесь в вики: https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Affine_transformation

Вам следует преобразовать F(X) в F(Y) следующим образом:

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)

Выходы:

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