Search code examples
pythonintegrationnormal-distribution

Python - Integrating entire list with the Gaussian distribution

I'm trying to plot the integral of this normal distribution for values greater than x defined as

enter image description here

where

enter image description here

defining both functions in python

import scipy.integrate as integrate
import numpy as np

def gaussian(x, mu, sig):
    norm = 1/np.sqrt(2*np.pi*sig*sig)
    return norm * np.exp(-np.power(x - mu, 2.) / (2. * sig*sig))

def gaussianGreater(x, mu, sig):
    Integrand = lambda x: gaussian(x, mu, sig)
    return integrate.quad(Integrand,-np.Inf, x)[0]

My problem now lies in the integration bounds of my gaussianGreater function while it is being evaluated through the distribution function. When evaluating, this occurs.

y = gaussianGreater(subdist_1, mu_1, sig_1 )
xd = np.argsort(subdist_1)

fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)
ax.plot(subdist_1[xd] ,y[xd] )

ValueError: The truth value of an array with more than one element
is ambiguous. Use a.any() or a.all()

I tried changing the upper bounds to what the error gave me, but that would return the error 'float' object has no attribute '__getitem__'

Applying a for loop does not work either

[gaussianGreater(x, mu_1, sig_1 ) for x in subdist_1]

TypeError: only integer arrays with one element can be converted to an index

How do I fix this problem?

Solution

  • You can directly use scipy.stats.norm's survival function for 1 - F(x):

    import scipy.stats as ss
x = np.linspace(-3, 3, 100)
y = ss.norm.sf(x)  # you can pass its mean and std. dev. as well
plt.plot(x, y)

    enter image description here