-2

I saw this formula: http://tutorial.math.lamar.edu/Classes/CalcI/ProofIntProp.aspx#Extras_IntPf_AvgVal and tried to implement python algorithm to approximate integrals. It kinda works, but does not make sense to me, so if any1 can explain why, it will be nice :) This is my code:

import random

def monte_carlo(function, a, b, iter = 100000):
    """
    function - 2d array of numbers, example: [[2, 1], [5, 4]] 2x^1 + 5x^4
    a, b - boundaries
    Approximates integral
    """ 
    answer = 0
    for i in range(0, iter): 
        rpt = random.randrange(a, b+1)
        print(i , 'th ' , 'iteration')
        answer += evall(function, rpt) 

    return (1/(b-a))*answer
def evall(function, point):
    result = 0
    for term in function:
        result += term[0]*pow(point, term[1])
    return result


print('Answer is: ', monte_carlo([[1, 2]], 1, 100))

and it works. But the formula says that we need the delta X in there, so if I make:

deltaX = (b-a)/iter

and then multiply evall(function, rpt) by it, it should work as well, but it does not. The example I used is for the function x^2.

Peter
  • 59
  • 1
  • 9
  • 2
    `...but it does not` - what is wrong with *it*? – wwii Feb 11 '20 at 22:35
  • Welcome to StackOverflow. [On topic](https://stackoverflow.com/help/on-topic), [how to ask](https://stackoverflow.com/help/how-to-ask), and ... [the perfect question](https://codeblog.jonskeet.uk/2010/08/29/writing-the-perfect-question/) apply here. StackOverflow is a knowledge base for *specific* programming problems -- not a tutorial resource. "Explain this code to me" is too broad. Show us your tracing of the code's operation, explain what you *do* understand, and the point where you get lost. – Prune Feb 11 '20 at 22:36
  • well I call it with [[1,2]] and boundaries 100 and 1, which is x^2 integrated between 1 and 100, which is 333 333, you can verify that the output would be far from that – Peter Feb 11 '20 at 22:36
  • shall I delete the post? – Peter Feb 11 '20 at 22:59

1 Answers1

1

Change your monte_carlo function to:

def monte_carlo(function, a, b, iter = 100000):
    """
    function - 2d array of numbers, example: [[2, 1], [5, 4]] 2x^1 + 5x^4
    a, b - boundaries
    Approximates integral
    """ 
    answer = 0
    for i in range(0, iter): 
        rpt = random.random()*(b-a) + a  # Change to continuous uniform
        print(i , 'th ' , 'iteration')   # Probably don't want this
        answer += evall(function, rpt) 

    return (b-a)*answer/iter             # Corrects the Monte Carlo part

Monte Carlo integration involves taking an average. You might have noticed that your approximation changes based on the number of iterations. The only thing that should change with the number of iterations is the Monte Carlo error.

mickey
  • 2,168
  • 2
  • 11
  • 20