scipy.integrate 错误:只有大小为 1 的数组可以转换为 Python 标量

问题描述

我试图通过改变标准偏差(std)来计算高斯分布的概率。 我期望通过使用具有 21 个点的高斯正交,在 -1 到 +1 的范围内积分,均值 = 0 并设置 std=1 和 std =2 将分别产生 p = 0.68 和 p = 0.95(附图片)。

import scipy.integrate as integrate
import math as m

#mean=0
#varying the sigma
def f(sigma,x):
    return m.exp(-1*(x**2)/(2*sigma**2))/(sigma*m.sqrt(2*m.pi))

def prob_at_nsigma(sigma):
    value = 0.
    
    ans,anserr =integrate.quadrature(f,-1,1,args=(sigma,),maxiter=21)
    value = ans
    return value

print(prob_at_nsigma(1))

并且我收到以下错误,我不明白为什么会出现“除以零”和“只有大小为 1 的数组可以转换为 Python 标量”:

runfile('C:/pythonExe/ntu_cp/untitled0.py',wdir='C:/pythonExe/ntu_cp')
C:\pythonExe\ntu_cp\untitled0.py:14: RuntimeWarning: divide by zero encountered in true_divide
  return m.exp(-1*(x**2)/(2*sigma**2))/(sigma*m.sqrt(2*m.pi))
C:\pythonExe\ntu_cp\untitled0.py:14: RuntimeWarning: invalid value encountered in true_divide
  return m.exp(-1*(x**2)/(2*sigma**2))/(sigma*m.sqrt(2*m.pi))
Traceback (most recent call last):

  File "C:\pythonExe\ntu_cp\untitled0.py",line 23,in <module>
    print(prob_at_nsigma(1))

  File "C:\pythonExe\ntu_cp\untitled0.py",line 19,in prob_at_nsigma
    ans,maxiter=21)

  File "C:\Users\cztee\anaconda3\lib\site-packages\scipy\integrate\_quadrature.py",line 238,in quadrature
    newval = fixed_quad(vfunc,a,b,(),n)[0]

  File "C:\Users\cztee\anaconda3\lib\site-packages\scipy\integrate\_quadrature.py",line 119,in fixed_quad
    return (b-a)/2.0 * np.sum(w*func(y,*args),axis=-1),None

  File "C:\Users\cztee\anaconda3\lib\site-packages\scipy\integrate\_quadrature.py",line 149,in vfunc
    return func(x,*args)

  File "C:\pythonExe\ntu_cp\untitled0.py",line 14,in f
    return m.exp(-1*(x**2)/(2*sigma**2))/(sigma*m.sqrt(2*m.pi))

TypeError: only size-1 arrays can be converted to Python scalars

感谢任何帮助。谢谢!

Gaussian distribution

解决方法

代码中有两个问题:

  • 如前所述,您应该使用 numpy 而不是 math 模块。
  • 您的 f 函数的参数顺序错误。 Scipy 假设第一个变量是 x 后跟其他参数。所以应该是f(x,sigma),而不是f(sigma,x)

解决方案:

import numpy as np
import scipy.integrate as integrate

def f(x,sigma):
    return np.exp(-1*(x**2)/(2*sigma**2))/(sigma*np.sqrt(2*np.pi))

def prob_at_nsigma(sigma):
    value = 0.
    ans,anserr = integrate.quadrature(f,-1,1,args=(sigma,),maxiter=21)
    return ans

print(prob_at_nsigma(1))
# 0.6826894922280757