使用 scipy.integrate.quad 时出现多个溢出警告

Multiple overflow warnings when using scipy.integrate.quad

我正在尝试在 Python 中实现以下功能:

r >= 0。代码是:

import scipy.integrate as integrate
from scipy.stats import norm
import numpy as np

def integrand(n, u, r):
    return norm.pdf(u, 0, 1) * ((norm.cdf(u + r, 0, 1) - 
                                 norm.cdf(u, 0, 1)) ** (n - 2)) * norm.pdf(u + r, 0, 1)

def fRn(n, r):
    return n * (n - 1) * integrate.quad(integrand, -np.inf, np.inf, args=(n, r))[0]

但是 fRn(2, 1) 的结果是:

__main__:14: RuntimeWarning: overflow encountered in double_scalars
__main__:17: IntegrationWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
Out[56]: nan

我找不到如何处理这个错误。我将不胜感激任何帮助。 谢谢


第二种方法

我试过以下方法:

from scipy.integrate import quad
import math
import numpy as np


def phi(x):
    return (1/math.sqrt(2*np.pi)) * np.exp(-x**2/2)

def PHI(x):
    return (1.0 + math.erf(x / math.sqrt(2.0))) / 2.0

def integrand(n, u, r):
    return phi(u) * ((PHI(u+r) - PHI(u))**(n-2)) * phi(u+r)


def fRn(n, r):
    return n * (n - 1) * quad(integrand, -np.inf, np.inf, args=(n, r))

而且,我收到以下错误:

Traceback (most recent call last):

  File "C:\Users\AppData\Local\Temp/ipykernel_15576/4053151620.py", line 1, in <module>
    fRn(2,0.1)

  File "C:\Users\Desktop\Simulation\Task_2.py", line 24, in fRn
    return quad(integrand, -np.inf, np.inf, args=(n, r))

  File "C:\Users\Anaconda3\lib\site-packages\scipy\integrate\quadpack.py", line 352, in quad
    points)

  File "C:\Users\Anaconda3\lib\site-packages\scipy\integrate\quadpack.py", line 465, in _quad
    return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

  File "C:\Users\Desktop\Simulation\Task_2.py", line 20, in integrand
    return phi(u) * ((PHI(u+r) - PHI(u))**(n-2)) * phi(u+r)

OverflowError: (34, 'Result too large')

积分变量应该是第一个参数,额外的参数(n, r)将在它之后传递,所以你的被积函数应该定义为

def integrand(u, n, r):
    return norm.pdf(u, 0, 1) * ((norm.cdf(u + r, 0, 1) - 
                                 norm.cdf(u, 0, 1)) ** (n - 2)
                             ) * norm.pdf(u + r, 0, 1)