How to create a mathematical expression in scipy?

Question

I have been exploring scipy and the core packages of scipy for a math project.

I need to do calculus operations on some equations... so for learning scipy I decided to test a simple equation (PDF of a normal random variable). I need to keep the constants as they are during the calculus operations... and not assign a value to it.

I was able to create it successfully using sympy. Here is the code:

from sympy import *

x = Symbol('x')
mu = Symbol('mu')
sigma = Symbol('sigma')

def normpdfeqn():
    y = exp(-(((x-mu)**2)) / (2*(sigma**2))) / (sqrt(2*pi*(sigma**2)))
    return y

print(integrate(normpdfeqn(), (x)))

and got a proper output:

sigma*erf(sqrt(2)*(-mu + x)/(2*sigma))/(2*sqrt(sigma**2))

Then I tried to do the same with scipy.

I have been reading http://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html but could not figure out how to create the equation for it. Here is what I have tried till now (it's almost the same as the above code):

from sympy import exp
from sympy import sqrt
from sympy import pi
from scipy.integrate import quad
from sympy import Symbol
x = Symbol('x')
mu = Symbol('mu')
sigma = Symbol('sigma')

def integrand():
    y = exp(-((x-mu)**2) / (2*(sigma**2))) / (sqrt(2*pi*(sigma**2)))
    return y


I = quad(integrand(), 0, 1,)
print(I)

The code might be far from complete, I do not know how to get it working.

If I would always be working with equations which I showed above, should I bother to learn scipy integrate? Or keep using sympy and numpy?

Answer 1

This simplified code might be a step towards your goals:

it takes a function, here a quadratic to keep things as simple as possible, and evaluates its closed form integral (using sympy)

Then it numerically evaluates the integral value between 0 and 3 using scipy.

And the closed form integral using the same limits with sympy.

It is a little bit of a prototype, but seems to do what you require.

from sympy import *
from scipy.integrate import quad

x = Symbol('x')

def func():
    y = x**2 + 2*x +3
    return y

def integrand(expr):
    y = lambda x: eval(expr)
    return y 


y = repr(func())             # a str representation of the original function

res = integrate(func(), (x)) # a sympy representation of the closed form integral 
                             # of the original function

I = quad(integrand(y), 0, 3) # a scipy numerical integration of the original function

print(y)        # sympy repr of original function
print(I)        # scipy numerical integration of the original function
print()
print(res)      # sympy repr of closed form integral of original function
print(integrand(repr(res))(3) - integrand(repr(res))(0)) # closed form integer evaluation

Resulting output:

x**2 + 2*x + 3
(27.0, 2.9976021664879227e-13)

x**3/3 + x**2 + 3*x
27.0