Numerical evaluation of a Gaussian Integral in Python?

Question

Goal

I'm trying to write code to compute the normalized Gaussian in the following,

$$ \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} $$

where $\mu \in [-10,10]$

Problem

Unfortunately, the integration algorithm does not converge and throws the warning:

FinalStat.py:68: IntegrationWarning: The integral is probably divergent, or slowly convergent.

The Gaussian should be normalizing to $\eqref{1}$, it seems like their an issue with Quadpack's backend. I'd like to have a fix for this to get the integral to normalize to its proper value?

Code

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)))

Output

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

Answer 1

In the example output from your code, $\sigma$ is huge, i.e. the Gaussian is extremely broad. The variable s you define as the pre-factor for the argument of the corresponding exponential is then only $\approx -1\cdot{}10^{-15}$, which is dangerously close to typical double precision limits (adding $10^{-16}$ to $1$ with typical double precision, e.g., still is $1$. scipy's quad having correspondingly to deal with huge and tiny numbers will make it hard to detect e.g. what numerically is zero.

For numerical integration to be stable, it is important to scale the integration variable appropriately: here, $\sigma$ is the typical length scale of your problem, and one would want the typical numerical scale used for integration to be of the order of $1$ (or $0.1$, or $10$, or anything of reasonable numerical magnitude). A possible substitution of the integration variable could be $u=(x-\mu)/\sigma$ (with ${\rm d}x=\sigma{\rm d}u$). The integral would then become: $$\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

Using the original integral over $f(x)$, I also find quad to fail and yield an incorrect integral of $\approx -4\cdot{}10^{-8}$ (using the same $\sigma$ as in the example output in the question, not another random value), employing $f2(u)$, quad yields $1$.

Let's say the transformation of the integrand shouldn't be as simple as above (which basically leads to a Gaussian with 'normalized' coefficients). One could also define $u=x/\sigma$, i.e. not absorb the shift due to $\mu$ in the transformation (which works here because $\mu \in [-10;10]$, if $\mu$ itself was huge, it would indeed be best to use the previous transformation to avoid numerical instabilities) and use the following integrand:

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

which, despite still having the shift due to $\mu$ in the integrand, accurately yields $1$ with scipy's quad. Choose e.g. $\mu=10^{12}$ and this latter approach will fail, quad again yielding incorrect results close to $0$.

A change of variables $u=x\cdot\alpha$ and ${\rm d}u={\rm d}x\cdot\alpha$ can generally be used for quadrature (and other numerical problems) to transform the integral kernel to numerically more tractable magnitudes.

Answer 2

No need for quadrature; your integral can be computed analytically. For this, you'll have to do the variable transformation $y = (x-\mu) / \sqrt{2} \sigma$: $$ \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 $$ (In your case, with random.seed(0), the result is -1.8286668122220665e-25 which explains the numerical difficulties you encouter.)