CODEWITHSUNDEEP
pythonintegrationnormal-distribution
iron2man
iron2man

Reputation: 1827

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?

Upvotes: 1

Views: 1427

Answers (2)

user2285236
user2285236

Reputation:

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

Upvotes: 1

spruceb
spruceb

Reputation: 621

What you've defined will work only when x is a single float. If what you want is to be able to pass a numpy array, say np.array([0.2, 0.3]) and get [N(>0.2), N(>0.3)], then you can use the function as you've already defined and just call the vectorize method. So if a was the numpy array you were using, a.vectorize(gaussianGreater) would get you an array with gaussianGreater applied to each element. You could define another function, vectorGausssianGreater that takes an array and just returns the result of that vectorize call on it. If what you really want is a single function that can take an array value, you could define something like the following:

def gaussian_greater(values, mean, standard_deviation):
    def greater_float(upper_bound):
        integrand = lambda y: gaussian(y, mean, standard_deviation)
        return integrate.quad(integrand, -np.inf, upper_bound)[0]
    return values.vectorize(greater)

Upvotes: 1

Related Questions