Python: Integrating a function but getting large errors

[MathewsMD]

My code:

import numpy as np
import matplotlib.pylab as plt
import math
import random
from scipy import integrate

R1 = .001
R2 = 7

def G(r,theta):
    sigma = random.randint(4000., 7000.)/1000. # width of beam is 4 - 7mm
    r0 = random.randint(0, R1*1000.)/1000. #random centroid
    theta0 = random.uniform(0, np.pi*2) #angle
    return (100/(np.pi*2*sigma**2))*(np.e**(-((r**2 + r0**2 - 2*r*r0*(np.cos(theta)*\
    np.cos(theta0) + np.sin(theta)*np.sin(theta0)))/2*sigma**2)))*r # this r is here because the integrand is dr and dtheta, NOT dx and dy any longer!
#amplitude is 100

RTI = integrate.nquad(G, [[0,R2],[0,2*np.pi]]) #Integration

print RTI

I keep running this code, and everything seems to be working, but I keep getting results like: (0.1406360958428128, 0.012357271987783056)

Also, the code itself takes quite a while to run and I'm sure this is an indication of something...just not sure where it's going on. As you can see, the error is almost the value of integration itself, and this error is to large. I've been looking at the code but am not sure why such a large error is arising. If anyone has any ideas, it would be greatly appreciated!

 

[phyzguy]

I'm not sure what exactly you are trying to do. The integrate.nquad function will call the function G many times as it evaluates the integral. Each time it calls G, the parameters sigma, r0 and theta0 will have different values. So the value of the integral isn't really defined, because the function G is a different function each time it is called. I suspect that this isn't what you intend. You should probably move the random.randint calls outside of the G function and call them once to define sigma, r0 and theta0 before evaluating the integral. Does this make sense?

 

[D H]

Also, the code itself takes quite a while to run and I'm sure this is an indication of something...

It's an indication that the integrator is struggling.

A lousy numerical integrator might arbitrarily chop the space to be integrated into fixed size chunks, call the derivative function at each grid point, and apply some mathematical formulae to generate an estimate of the integral over than space. A better integrator will subdivide the space according to sensitivities revealed by calls to the derivative function.

Your derivative function involves random. A good integrator will detect that. What a good integrator won't do is detect that your derivative function (in the context of computer science) is not a "function" (in the context of mathematics). Any integrator must necessarily assume that the function to be integrated is indeed a function. You don't have a function. You have randomness.

As you can see, the error is almost the value of integration itself, and this error is to large.

What, exactly, are you trying to integrate?

 

[MathewsMD]

I'm not sure what exactly you are trying to do. The integrate.nquad function will call the function G many times as it evaluates the integral. Each time it calls G, the parameters sigma, r0 and theta0 will have different values. So the value of the integral isn't really defined, because the function G is a different function each time it is called. I suspect that this isn't what you intend. You should probably move the random.randint calls outside of the G function and call them once to define sigma, r0 and theta0 before evaluating the integral. Does this make sense?

I have attached an image of what I am trying to do. Yes, you seem to get what I was attempting. I want to generate a single particular gaussian with the variables r0, theta0 and sigma0 randomly generated but held constant for that one run, then I want to integrate that one particular function over the circle. I would be creating numerous random gaussian functions like this.

I have modified the code a little bit based on your suggestions, and it seems to be working, but any input with regards to thingsto modify would be great!

Code:

RTIT = []

R1 = 4. #
R2 = 7. #

i = 0

    while i < 10:

    i = i + 1

    sigma0 = random.randint(4000., 7000.)/1000. # width is 4 - 7mm
    r0 = random.randint(0, R1*1000.)/1000. #random centroid
    theta0 = random.uniform(0, np.pi*2) #angle centroid makes with polar axis

    def G(r,theta): #this is the Gaussian
        return (100/(np.pi*2*sigma0**2))*(np.e**(-((r**2 + r0**2 \
        - 2*r*r0*(np.cos(theta)*np.cos(theta0) + \
        np.sin(theta)*np.sin(theta0)))/2*sigma0**2)))*r # this r is here because
#the integrand is dr and dtheta, NOT dx and dy any longer!
#intensity is 100

    RTI = integrate.nquad(G, [[0,R2],[0,2*np.pi]]) #Total
    RTIT.append(RTI)

print RTIT

 

[phyzguy]

A couple of things:

If you enclose your code in

tags, then the code will be more readable.

Your code with the modifications looks good. Did it fix the large error problem you were initially having?

 

[MathewsMD]

A couple of things:

If you enclose your code in

tags, then the code will be more readable.

Your code with the modifications looks good. Did it fix the large error problem you were initially having?

Yep, it seems to be working great now. Thanks! I'll definitely use that from now on.