目标

我正在尝试编写代码来计算下面的归一化高斯，

\begin{matrix} (1) & \int_{- \infty}^{\infty} \frac{1}{σ \sqrt{2 π}} \exp (- \frac{(x - μ)^{2}}{2 σ^{2}}) d x \end{matrix}

$\begin{equation} \int_{-\infty}^{\infty} \frac{1}{ \sigma \sqrt{2 \pi}} \exp\bigg( - \frac{(x - \mu)^{2}}{2 \sigma^{2}}\bigg)dx \label{1} \tag{1} \end{equation}$

在哪里 $\mu \in [-10,10]$

问题

不幸的是，积分算法不收敛并抛出警告：

FinalStat.py:68：IntegrationWarning：积分可能是发散的，或缓慢收敛的。

高斯应该归一化为 $\eqref{1}$ ，这似乎是 Quadpack 后端的问题。我想对此进行修复以使积分正常化为其正确值？

代码

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

xvalues = []
yvalues = []
def generate():
    #=================================================================== #
    #                                                                    #
    #                 Generates Linear Data                              #
    #                a,b are random varibles                             #
    # ================================================================== #
    for i in range(0,10):
        a = random.randint(-10, 10)
        b = random.randint(-10, 10)
        xvalues.append(i)
        y = a * (b + i)
        yvalues.append(y)


def weighted_mean(yvalues):
    #=============================================#
    #          Computes the Weighted Mean         #
    #=============================================#
    y_i = np.array(yvalues)
    x_i = np.array(xvalues)
    evaulated_mean = sum(x_i*y_i) / len(y_i)
    return evaulated_mean

def weighted_variance(yvalues):
    #============================================#
    #       Computes the Weighted Variance       #
    #============================================#

    s =  []

    yvalues_mean = weighted_mean(yvalues)

    for y in yvalues:
        z = (y - yvalues_mean)**2
        s.append(z)
    #===================================#
    #s_i = value of the data set        #
    #v_i = Number of Data Points in Pop #
    #===================================#
    s_i = np.array(s)
    v_i = np.array(s)
    t = (sum(s_i*v_i)/len(v_i)) 
    print("The Variance", t)
    return t



def gaussian(sigma,mu, x):
    #===================================================#
    #Define and Compute Gaussian Function with the FWDM #
    #===================================================#   
    FWHM = 2*(numpy.sqrt(2*numpy.log(2)))*sigma
    k = 1 / (sigma * math.sqrt(2*math.pi))
    s = -1.0 / (2 * sigma * sigma)
    def f(x):
        return k * math.exp(s * (x - mu)*(x - mu))

    print("The corresponding FHWM", FWHM)
    print("The Integral is", quad(f, -np.inf, np.inf))
    return FWHM


generate()
print( "The Mean is =>" , weighted_mean(yvalues))
#weighted_variance(yvalues)
print("our Normal Distrubtion Equals" , gaussian(weighted_variance(yvalues), weighted_mean(yvalues), random.randint(-10,10)))

输出

The Mean is => -53.7
The Variance 18538654.5401
The corresponding FHWM 43655195.3189315
FinalStat.py:68: IntegrationWarning: The integral is probably divergent, or slowly convergent.print("The Integral is", quad(f, -np.inf, np.inf))
The Integral is (-4.3038967971581716e-08, 6.997796078221529e-13)
our Normal Distrubtion Equals 43655195.3189315

2个回答

无需正交；您的积分可以解析计算。为此，您必须进行变量转换 $y = (x-\mu) / \sqrt{2} \sigma$ ：

\int_{- \infty}^{\infty} k \exp (- \frac{(x - μ)^{2}}{2 σ^{2}}) d x = k \sqrt{2} σ \int_{- \infty}^{\infty} \exp (- y^{2}) d y = \sqrt{2 π} k σ

$\int\limits_{-\infty}^{\infty} k\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) {\rm d}x = k\sqrt{2}\sigma \int\limits_{-\infty}^{\infty} \exp\left(-y^2\right) {\rm d}y = \sqrt{2\pi} k\sigma$ （在您的情况下，random.seed(0)结果是-1.8286668122220665e-25解释了您遇到的数字困难。）

在代码的示例输出中， $\sigma$ 是巨大的，即高斯非常广泛。s您定义为相应指数参数的前置因子的变量只有 $\approx -1\cdot{}10^{-15}$ ，这危险地接近典型的双精度限制（添加 $10^{-16}$ 到 $1$ 具有典型的双精度，例如，仍然是 $1$ . scipy's quad相应地处理巨大和微小的数字将使得难以检测例如数字为零。

为了使数值积分稳定，适当缩放积分变量很重要：这里， $\sigma$ 是您的问题的典型长度尺度，并且希望用于积分的典型数值尺度为 $1$ （或者 $0.1$ ，要么 $10$ ，或任何合理的数值大小）。积分变量的可能替换可能是 $u=(x-\mu)/\sigma$ （和 ${\rm d}x=\sigma{\rm d}u$ ）。然后积分将变为：

\int_{- \infty}^{\infty} k \exp (- \frac{(x - μ)^{2}}{2 σ^{2}}) d x = \int_{- \infty}^{\infty} σ k \exp (- \frac{u^{2}}{2}) d u .

$\int\limits_{-\infty}^{\infty} k\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) {\rm d}x = \int\limits_{-\infty}^{\infty} \sigma k\exp\left(-\frac{u^2}{2}\right) {\rm d}u .$

def f2(u):
    return k*numpy.exp(-0.5*u**2)*sigma

使用原始积分 $f(x)$ ，我也发现quad失败并产生不正确的积分 $\approx -4\cdot{}10^{-8}$ （使用相同的 $\sigma$ 如问题中的示例输出所示，而不是另一个随机值），采用 $f2(u)$ ,quad产量 $1$ .

假设被积函数的变换不应该像上面那样简单（这基本上会导致具有“归一化”系数的高斯）。也可以定义 $u=x/\sigma$ ，即不吸收由于 $\mu$ 在转换中（在这里有效，因为 $\mu \in [-10;10]$ ，如果 $\mu$ 本身很大，确实最好使用前面的转换来避免数值不稳定）并使用以下被积函数：

def f3(u):
   return k*numpy.exp(-0.5*(u*sigma-mu)**2/sigma**2)*sigma

其中，尽管由于 $\mu$ 在被积函数中，准确地产生 $1$ 与scipy's quad. 选择例如 $\mu=10^{12}$ 后一种方法将失败，quad再次产生接近于 $0$ .

变量的变化 $u=x\cdot\alpha$ 和 ${\rm d}u={\rm d}x\cdot\alpha$ 通常可用于求积（和其他数值问题）以将积分核转换为数值上更易于处理的量级。

其它你可能感兴趣的问题

上一篇当我使用具有对称半正定矩阵的共轭梯度求解器时会发生什么？下一篇为什么矫顽力在 FEM 框架中如此重要？