2

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.

Ignacio Vergara Kausel
  • 5,521
  • 4
  • 31
  • 41
Curious
  • 33
  • 1
  • 7
  • Concerning `sympy` you should help it by properly defining the type of the symbols. Most likely n is an integer/natural while z should be positive real. As for the numerical integration, `scipy.stats.chi2` returns an object and not a value. You have to find out how to get a value of the `chi2` distribution evaluated at the z you want. – Ignacio Vergara Kausel Feb 02 '18 at 11:16
  • For the numerical integration, i do not know how to evaluate chi2 of z because z is a symbol... N is a natural, while z is complex I would say. But even when I try to define n as natural and z as integer (although i think it is complex), Kernel still runs for eternity... – Curious Feb 02 '18 at 12:18
  • z is a number, which for your D function is the integration variable. Without giving `sympy` hints, helps, and splitting the problem in pieces it'll have hard time figuring out a solution. You should really try for numerical integration, worst case you program your chi2 function. – Ignacio Vergara Kausel Feb 02 '18 at 12:43

2 Answers2

3

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

1

For integrating equations having symbols other than the variable of integration you have to use sympy library.

What you have to do is actually define the symbols as symbol (Python doesn't knows if those alphabets(symbols) represent some real valued numbers and are meant to be left out during integration).

After integration is done you may substitute the values of symbols with some numbers.

Please check the link below for example code of integrating equations with symbols

EXAMPLE CODE

Check the answer by user : MSD in the link, it will make it clear

MSD
  • 154
  • 2
  • 12