Python Bivariate Normal CDF с переменной верхней границей

Question

Я пытаюсь найти элегантный способ вычисления двумерного нормального CDF с питоном, где одна верхняя граница CDF является функцией двух переменных, одной из которых является переменная двумерной нормальная плотность (интегральная величина).

Пример:

from scipy import integrate 
import numpy as np 

# First define f(x, y) as the bivariate normal distribution with fixed correlation p 
p = 0.4 
def f(x, y): 
    Q = x**2 + y**2 - 2*p*x*y 
    return 1/(2*np.pi*np.sqrt(1-p**2))*np.exp(-1/(2*(1-p**2))*Q) 

# Define N2(a, b) as the cumulative bivariate normal distribution f where a is constant 
# and b is a function of two variables 
def N2(a, b): 
    prob, error = integrate.dblquad(f, np.NINF, a, lambda x: np.NINF, b) 
    return prob 

# Upper bound function of 2 variables example where x is an integral variable 
def upper_bound(x, v): 
    return 0.5*v*x 

# My approach which doesn't work 
# Calculate bivariate normal CDF for different values of v in the upper bound 
results = [N2(1, upper_bound(x, v)) for v in np.linspace(0.01, 4, 200)] 

Любые идеи о том, как я мог бы изменить свой подход, поэтому призыв к upper_bound(x, v) в results = [N2(1, upper_bound(x, v)) for v in np.linspace(0.01, 4, 200)] будет работать? Также приветствуются другие подходы к решению этой проблемы.

Редактировать: Это интеграл, который я хочу вычислить, где f (x, y) - бинарная нормальная функция плотности. Обратите внимание, что фактическая верхняя граница f (x, v) = 0,5 * v * x, которую я хочу вычислить, является более сложной, это как раз пример, поэтому я не хочу его вычислять символически, например, с помощью sympy. Кроме того, моя цель состоит в том, чтобы вычислить интеграл для нескольких сотен различных значений V

Интеграл:.

Answer 1

Хотя это медленный такой подход, кажется, работает.

Первые несколько строк, до «это должно произвести 1», являются проверкой работоспособности. Я хотел убедиться, что мой подход будет правильно рассчитать объем под плотностью. Оно делает.

Я использую матрицу дисперсии-ковариации, чтобы получить требуемую корреляцию 0,4 и не писать собственный PDF-файл.

Я выполняю функции в двух местах, чтобы функции имели только одиночные параметры. Это позволяет вычислить внутренний интеграл как функцию x. Это также позволяет взять параметр «вне» других расчетов.

from toolz import curry 
from scipy.stats import multivariate_normal 
from scipy.integrate import quad 
import numpy as np 

@curry 
def bivariate(x,y): 
    return multivariate_normal.pdf([x,y],cov=[[1,.4],[.4,1]]) 

def out_y(x): 
    marginal = bivariate(x) 
    return quad(marginal, np.NINF, np.PINF)[0] 

# this should produce 1 
print (quad(out_y, np.NINF, np.PINF)[0]) 

# now to what the OP wants 

@curry 
def inner_integral(v,x): 
    marginal = bivariate(x) 
    return quad(marginal, np.NINF, 0.5*v*x)[0] 

inner_integral_for_one_v = inner_integral(0.8) 
print (quad(inner_integral_for_one_v,np.NINF, 1)[0]) 

Чтобы использовать этот код вы бы написать что-то эквивалентно:

for v in range(0,1,0.1): 
    inner_integral_for_one_v = inner_integral(v) 
    print (quad(inner_integral_for_one_v,np.NINF, 1)[0])