Search code examples
pythonscipysympynumerical-integration

Numerical integration of symbolic function in Python

I am new to Python so some of my questions or ideas might be silly, but...

I want to graph a distribution D(x). m and s2 are some given real numbers. I have been told that the best way to graph D(x) is to write a function which solves the integral (which is inside a function D(x)) for every x.

D(x) function

where chi2 is so defined:

enter image description here

So, as much as I know math, I am supposed to integrate first and then I can solve for each x, correct me if I am wrong.

I have also been told to calculate the integral numerically, but I do not know how to do it because the functions includes symbols. I have already tried using symbolical integration (despite what I have been told), but the kernel never ends the process of integration, as well as when I try to compute it numerically. When I tried integrationg numericaly, I used lamdify, of course.

So here are my codes: 1. trying to solve symbolic integral

from sympy import symbols, integrate, sqrt, exp, oo

s2= 0.0628777415586
m= 5.02422436191

x, n, z=symbols ('x, n, z')
integrate(exp(-n/(z+1) * (x-m)**2/2*s2)  *  1/2 / sqrt(z+1), (z, 0, oo))

Not working, Kernel never stops. [I put 1/2 instead of chi2 in formula, intented to change it later]

An alternative is trying to solve numerical integral (calling it from a function D(x))

import scipy.integrate
from scipy.stats import chi2
from math import *


s2= 0.0628777415586
m= 5.02422436191

def numint(z, n, x):
    return exp(-n/(z+1) * (x-m)**2/2*s2)  *  scipy.stats.chi2(z) / sqrt(z+1)

scipy.integrate.quad(numint, 0, np.inf, args=(1, 2))

Error comes from the line with return and says: 

TypeError: unsupported operand type(s) for *: 'float' and 'rv_frozen'

I suppose that the problem here comes becaue of chi2 which is rv_frozen, but how to make it work? Is there anything in my code that is wrong? I have no idea if it is right what I wrote and how to fix this... I have been working on this for a very long time and am a bit desperate, so any help is welcome.

Solution

  • I am supposed to integrate first and then I can solve for each x, correct me if I am wrong.

    No, you pick x and then integrate. The only chance to integrate this function is numerically, and for that the values of all symbols except the variable of integration (z) must be given.

    m and s2 are some given real numbers.

    Then give them. Also, n must be specified. The function D(x) depends on symbolic parameters m, n, s (not z, because z is being integrated out). Graphing requires numeric evaluation, and this function cannot be evaluated numerically until some values of m, n, s are chosen.

    the problem here comes becaue of chi2 which is rv_frozen

    Yes, chi2 is a random variable and you want the probability density function of that random variable, accessed by .pdf method. Read the docs.

    A complete example.

    import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.stats import chi2
integrand = lambda z, x, m, n, s: np.exp(-(n/(z+1))*(x-m)**2/(2*s**2)) * chi2.pdf(z, n)/np.sqrt(z+1)
D = lambda x, m, n, s: np.sqrt(n/2*np.pi*s**2) * quad(integrand, 0, np.inf, args=(x, m, n, s))[0]
x = np.linspace(0, 5, 100)
y = np.vectorize(D)(x, 3, 10, 1.5)  # some values of m, n, s
plt.plot(x, y)
plt.show()

    graph