5
$\begingroup$

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
$\endgroup$
  • 1
    $\begingroup$ I think trapezoidal quadrature rather than the default from scipy.integrate.quad will help. $\endgroup$ – user14717 Jun 15 at 16:22
  • 3
    $\begingroup$ Might want to have a look at Gauss-Hermite quadrature: en.wikipedia.org/wiki/Gauss%E2%80%93Hermite_quadrature $\endgroup$ – sssssssssssss Jun 15 at 16:35
  • 2
    $\begingroup$ Wouldn't all the integrals give the same value that, by the way, can be computed analytically. Or, am I missing something? $\endgroup$ – nicoguaro Jun 15 at 18:18
  • 2
    $\begingroup$ @nicoguaro: I figured this was just an exercise; not a claim that this is the best way to do it. $\endgroup$ – user14717 Jun 15 at 19:40
  • 2
    $\begingroup$ In that case, I would also suggest to use Hermite quadrature. $\endgroup$ – nicoguaro Jun 15 at 19:52
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I didn't know that SciPy's quad had such limits also could turning the Gaussian into an IVP problem help cure the issue $\endgroup$ – Zophikel Jun 17 at 13:45
4
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Could you give more details about the random.seed(0), and how it relates to the result ? $\endgroup$ – Zophikel Jun 17 at 14:56
  • $\begingroup$ The computation uses randomness (in generate()), so whatever the values are there, they will affect the final result. Fixing the random.seed makes the results reproducible. $\endgroup$ – Nico Schlömer Jun 17 at 17:21
  • 1
    $\begingroup$ One thing is that the script will yield different results every time since random.seed is not fixed. The other is that one can test whether scipy.integrate.quad can handle Gaussian input with any randomly chosen width $\sigma$: the function definition with the choice for the pre-factor $k$ indeed uses the fact that there is an analytical result of the integral for normalization. The automatic integrator quad fails, if $\sigma$ is too large. The numerical difficulties aren't really the randomness of the input itself, but rather that a random width can be too large for quad to handle. $\endgroup$ – v-joe Jun 17 at 18:46
  • 2
    $\begingroup$ @v-joe Not really. The variance depends on the seed, but it is always too large. $\endgroup$ – Nico Schlömer Jun 17 at 20:35
  • 1
    $\begingroup$ @NicoSchlömer Good point. Looking at it that way (i.e. generating random data and computing the variance of that) there is really no need to employ any notion of a Gaussian = (continuous) normal distribution and any quadrature whatsoever but just take the variance computed using a discrete sum. $\endgroup$ – v-joe Jun 17 at 20:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.