CODEWITHSUNDEEP
pythonarraysscipynormal-distributionvariance
Duc Anh Ho
Duc Anh Ho

Reputation: 21

how to integrate gaussian distribution with an array of variance

I'm trying to calculate a integration of gaussian distribution. I have an array of sigma.

sigma = np.array([0.2549833 , 0., 0.42156247, 0. , 0.,  0., 0.79124217, 0.54235005, 0.79124217, 0.        ,     0.        , 0.32532629, 0.46753655, 0.60605513, 0.55420338,      0.        , 0.38053264, 0.42690288, 0.        , 0.63526099])

And the formula of gaussian distribution:

def gaussian(x, mu, sig):
if sig != 0:
    return np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))

Integrate this gaussian distribution:

I = np.zeros(len(sigma), dtype=float)
for i in range(0, len(sigma)):
    I[i] = quad(gaussian(x, mu = 0, sig = sigma[i]), 0, 105)

but it's not working because quad function is giving error. How can I get an array of Integration in this case?

Upvotes: 0

Views: 385

Answers (2)

ev-br
ev-br

Reputation: 26090

Integration of the Gaussian pdf gives the cumulative density function, cdf. So, what you need is scipy.stats.norm.cdf. And that one supports vectorized evaluations, so no need for a for-loop in your code.

Note that cdf integrates from the negative infinity. Also, the sigma parameter is called scale.

Upvotes: 0

Sheldore
Sheldore

Reputation: 39072

The problem was in your usage of quad. Here is the correct version. Few things to note:

  • You are creating I of length equal to number of sigma values, but you were doing nothing when sigma is 0. So now, I recreated sigma to have only non-zero values. This way I removed the condition if sig != 0:

  • The correct usage as per docs is that you first pass the function definition, then the limits and then are function arguments using the keyword args.

  • The quad will return a tuple containing the integral and the absolute error. Hence, you should use [0] index to get the integral value to store in I. You can store the absolute error in a different list if you want.

I added the plot of Integral versus sigma values for your convenience.

Modified part of your code

sigma = sigma[sigma !=0]

def gaussian(x, mu, sig):
        return np.exp(-(x - mu)**2/ (2 * sig**2))

I = np.zeros(len(sigma), dtype=float)

for i in range(0, len(sigma)):
    I[i]  = quad(gaussian, 0, 105,  args=(0, sigma[i]))[0]

enter image description here

Upvotes: 1

Related Questions