Python 具有可变上限的双变量正态 CDF
Python Bivariate Normal CDF with variable upper bound
我正在尝试找到一种优雅的方法来计算具有 python 的双变量正态 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 值的积分。
积分：
虽然速度很慢，但这种方法似乎有效。
前几行，直到 'this should produce 1'，是完整性检查。我想验证我的方法是否能正确计算密度下的体积。确实如此。
我使用 variance-covariance 矩阵来获得所需的 0.4 相关性，避免编写自己的 pdf。
我在两个地方柯里化函数，这样函数只有一个参数。这使得可以将内积分计算为 x 的函数。也使得取v参数'outside'其他计算成为可能。
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])
我不得不在 Python 中编写一个使用双变量分布的期权模型。但是，我没有找到一个快速的预建函数 - 有些似乎使用随机 scipy 生成器来模拟它与多元函数。但是......如果你真的深入挖掘并查看其他财务软件包正在使用什么，那么它几乎都是由华盛顿大学的 Alan Genz 一个人编写的所有代码。他几乎用 Fortan 或 MATLAB 编写所有内容，仅此而已。因此，如果您查看具有双变量 CDF 的包，您会在那里找到他的名字和他的代码（实际上我在 MATLAB 中找到了一个旧版本）。 http://www.math.wsu.edu/faculty/genz/software/software.html
那么为什么 SciPy 或 NumPy 还没有内置它，我不知道。但我在几个小时内重写了它，同时使用 MATLAB 和 Python 检查生成的代码。他正在用高斯勒让德正交求解，所以这总是比使用随机数生成器的解决方案快得多：https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature
# Alan Genz - original MATLAB code converted by Matt Slezak to Python
# http://www.math.wsu.edu/faculty/genz/software/matlab/bvnl.m
# Department of Mathematics
# Washington State University
# This is the bivariate CDF distribution
#
# dh 1st upper integration limit
# dk 2nd upper integration limit
# r correlation coefficient
import numpy as np
from scipy.stats import norm
def bivariate_cdf(dh, dk, r):
# dh and dk get signs flipped right away
dh = float(-dh)
dk = float(-dk)
if (dh == np.inf) or (dk == np.inf):
return 0
else:
if dh == - np.inf:
if dk == - np.inf:
return 1
else:
return norm.cdf(- dk)
if dk == - np.inf:
return norm.cdf(- dh)
else:
if r == 0:
return norm.cdf(- dh)*norm.cdf(- dk)
else:
tp=2*np.pi
h=dh
k=dk
hk=h*k
bvn=0
if abs(r) < 0.3:
w=np.array([0.1713244923791705,0.3607615730481384,0.4679139345726904])
x=np.array([0.9324695142031522,0.6612093864662647,0.238619186083197])
else:
if abs(r) < 0.75:
w=np.array([0.04717533638651177,0.1069393259953183,0.1600783285433464,0.2031674267230659,0.2334925365383547,0.2491470458134029])
x=np.array([0.9815606342467191,0.904117256370475,0.769902674194305,0.5873179542866171,0.3678314989981802,0.1252334085114692])
else:
w=np.array([0.01761400713915212,0.04060142980038694,0.06267204833410905,0.08327674157670475,0.1019301198172404,0.1181945319615184,0.1316886384491766,0.1420961093183821,0.1491729864726037,0.1527533871307259])
x=np.array([0.9931285991850949, 0.9639719272779138, 0.9122344282513259,0.8391169718222188, 0.7463319064601508, 0.6360536807265150,0.5108670019508271,0.3737060887154196,0.2277858511416451,0.07652652113349733])
w = np.concatenate((w, w), axis=0)
x = np.concatenate((1 - x, 1 + x), axis=0)
if abs(r) < 0.925:
hs=(h*h+k*k) / 2
asr=np.arcsin(r) / 2
sn=np.sin(asr*x)
bvn=np.dot(np.exp((sn*hk-hs)/(1-sn**2)),w.T)
bvn=bvn*asr/tp + norm.cdf(-h)*norm.cdf(-k)
else:
if r < 0:
k=- k
hk=- hk
if abs(r) < 1:
as1=1 - r ** 2
a=np.sqrt(as1)
bs=(h - k) ** 2
asr=- (bs / as1 + hk) / 2
c=(4 - hk) / 8
d=(12 - hk) / 80
if asr > - 100:
bvn= a*np.exp(asr)*(1-c*(bs-as1)*(1-d*bs)/3+c*d*as1**2)
if hk > - 100:
b=np.sqrt(bs)
sp=np.sqrt(tp)*norm.cdf(-b/a)
bvn=bvn - np.exp(-hk/2)*sp*b*( 1 - c*bs*(1-d*bs)/3 )
a=a / 2
xs=(a*x) ** 2
asr=- (bs / xs + hk) / 2
ix=asr > - 100
xs=xs[ix]
sp=( 1 + c*xs*(1+5*d*xs) )
rs=np.sqrt(1 - xs)
ep=np.exp(- (hk / 2)*xs / (1 + rs) ** 2) / rs
bvn=( a*( np.dot((np.exp(asr[ix])*(sp-ep)),w[ix].T) ) - bvn )/tp
if r > 0:
bvn=bvn + norm.cdf(- max(h,k))
else:
if h >= k:
bvn=- bvn
else:
if h < 0:
L=norm.cdf(k) - norm.cdf(h)
else:
L=norm.cdf(- h) - norm.cdf(- k)
bvn=L - bvn
return max(0,min(1,bvn))
我正在尝试找到一种优雅的方法来计算具有 python 的双变量正态 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 值的积分。
积分：
虽然速度很慢，但这种方法似乎有效。
前几行，直到 'this should produce 1'，是完整性检查。我想验证我的方法是否能正确计算密度下的体积。确实如此。
我使用 variance-covariance 矩阵来获得所需的 0.4 相关性，避免编写自己的 pdf。
我在两个地方柯里化函数，这样函数只有一个参数。这使得可以将内积分计算为 x 的函数。也使得取v参数'outside'其他计算成为可能。
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])
我不得不在 Python 中编写一个使用双变量分布的期权模型。但是，我没有找到一个快速的预建函数 - 有些似乎使用随机 scipy 生成器来模拟它与多元函数。但是......如果你真的深入挖掘并查看其他财务软件包正在使用什么，那么它几乎都是由华盛顿大学的 Alan Genz 一个人编写的所有代码。他几乎用 Fortan 或 MATLAB 编写所有内容，仅此而已。因此，如果您查看具有双变量 CDF 的包，您会在那里找到他的名字和他的代码（实际上我在 MATLAB 中找到了一个旧版本）。 http://www.math.wsu.edu/faculty/genz/software/software.html
那么为什么 SciPy 或 NumPy 还没有内置它，我不知道。但我在几个小时内重写了它，同时使用 MATLAB 和 Python 检查生成的代码。他正在用高斯勒让德正交求解，所以这总是比使用随机数生成器的解决方案快得多：https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature
# Alan Genz - original MATLAB code converted by Matt Slezak to Python
# http://www.math.wsu.edu/faculty/genz/software/matlab/bvnl.m
# Department of Mathematics
# Washington State University
# This is the bivariate CDF distribution
#
# dh 1st upper integration limit
# dk 2nd upper integration limit
# r correlation coefficient
import numpy as np
from scipy.stats import norm
def bivariate_cdf(dh, dk, r):
# dh and dk get signs flipped right away
dh = float(-dh)
dk = float(-dk)
if (dh == np.inf) or (dk == np.inf):
return 0
else:
if dh == - np.inf:
if dk == - np.inf:
return 1
else:
return norm.cdf(- dk)
if dk == - np.inf:
return norm.cdf(- dh)
else:
if r == 0:
return norm.cdf(- dh)*norm.cdf(- dk)
else:
tp=2*np.pi
h=dh
k=dk
hk=h*k
bvn=0
if abs(r) < 0.3:
w=np.array([0.1713244923791705,0.3607615730481384,0.4679139345726904])
x=np.array([0.9324695142031522,0.6612093864662647,0.238619186083197])
else:
if abs(r) < 0.75:
w=np.array([0.04717533638651177,0.1069393259953183,0.1600783285433464,0.2031674267230659,0.2334925365383547,0.2491470458134029])
x=np.array([0.9815606342467191,0.904117256370475,0.769902674194305,0.5873179542866171,0.3678314989981802,0.1252334085114692])
else:
w=np.array([0.01761400713915212,0.04060142980038694,0.06267204833410905,0.08327674157670475,0.1019301198172404,0.1181945319615184,0.1316886384491766,0.1420961093183821,0.1491729864726037,0.1527533871307259])
x=np.array([0.9931285991850949, 0.9639719272779138, 0.9122344282513259,0.8391169718222188, 0.7463319064601508, 0.6360536807265150,0.5108670019508271,0.3737060887154196,0.2277858511416451,0.07652652113349733])
w = np.concatenate((w, w), axis=0)
x = np.concatenate((1 - x, 1 + x), axis=0)
if abs(r) < 0.925:
hs=(h*h+k*k) / 2
asr=np.arcsin(r) / 2
sn=np.sin(asr*x)
bvn=np.dot(np.exp((sn*hk-hs)/(1-sn**2)),w.T)
bvn=bvn*asr/tp + norm.cdf(-h)*norm.cdf(-k)
else:
if r < 0:
k=- k
hk=- hk
if abs(r) < 1:
as1=1 - r ** 2
a=np.sqrt(as1)
bs=(h - k) ** 2
asr=- (bs / as1 + hk) / 2
c=(4 - hk) / 8
d=(12 - hk) / 80
if asr > - 100:
bvn= a*np.exp(asr)*(1-c*(bs-as1)*(1-d*bs)/3+c*d*as1**2)
if hk > - 100:
b=np.sqrt(bs)
sp=np.sqrt(tp)*norm.cdf(-b/a)
bvn=bvn - np.exp(-hk/2)*sp*b*( 1 - c*bs*(1-d*bs)/3 )
a=a / 2
xs=(a*x) ** 2
asr=- (bs / xs + hk) / 2
ix=asr > - 100
xs=xs[ix]
sp=( 1 + c*xs*(1+5*d*xs) )
rs=np.sqrt(1 - xs)
ep=np.exp(- (hk / 2)*xs / (1 + rs) ** 2) / rs
bvn=( a*( np.dot((np.exp(asr[ix])*(sp-ep)),w[ix].T) ) - bvn )/tp
if r > 0:
bvn=bvn + norm.cdf(- max(h,k))
else:
if h >= k:
bvn=- bvn
else:
if h < 0:
L=norm.cdf(k) - norm.cdf(h)
else:
L=norm.cdf(- h) - norm.cdf(- k)
bvn=L - bvn
return max(0,min(1,bvn))