CODEWITHSUNDEEP
pythonscipystatisticsnumerical-methodsnumerical-integration
Sanit
Sanit

Reputation: 80

Calculating expectation of functions across normal distribution

I want to compute the expectation of certain functions across a normal distribution.

An example:

mu = 100
k = 100
sigma = 10
val, err = quad(lambda x: norm.pdf((x - mu) / sigma) * x if x > k else 0, -math.inf, math.inf)
print(val)

This prints 4.878683842492743e-288 which is clearly not the correct answer.

I assume this is happening because SciPy is unable to integrate the Gaussian. How can I solve this? Ideally, I'd want a method that'd allows one to integrate all sorts of functions across Gaussian and is not specific to the example I have put in.

Thanks!

Upvotes: 1

Views: 867

Answers (1)

dkaramit
dkaramit

Reputation: 125

I think this is a problem with the quadrature (sometimes it doesn't really adapt), and doesn't like the if statement.

So I would suggest something like this (integrate from k to infinity):

def f(x):
    return 1/sigma*norm.pdf((x - mu) / sigma)*x

val, err = quad(f,  k, math.inf)

Notice, as implied by Jimmy, the correct form of the Gaussian needs 1/sigma.

Another way to do this integral would be to force quad to be careful at some points. My favorite way is to do something like

import numpy as np
from scipy.integrate import quad 

#this is the Gaussian. Note that *0.5*(np.sign(x-k)+1) is 0 for x<k and 1 otherwise.
def f(x):
    return 1/(np.sqrt(2*np.pi)*sigma)*np.exp(- 0.5*( (x-mu)/sigma )**2 ) *x*0.5*(np.sign(x-k)+1)
    
#use this to integrate from in (-1,1)
def G(u):
    x=u/(1-u**2)
    return f(x)*(1+u**2)/(1-u**2)**2

quad(G,-1,1,points=np.linspace(-0.999,0.999,25))

I would suggest to read this in order to get how you can optimize such integrals.

Upvotes: 1

Related Questions