<jupyter_start><jupyter_text>Monte Carlo Integration by <NAME> ASTR 119  Here I will delcare all the necessary data packs for the following code to work. Other than the usual Numpy and Matplotlib, I will be using scipy's 'random' library for its built in random number generator.<jupyter_code>from scipy import random
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Here is our function of interest cos(x). We will keep this as a function so we can repeatedly call it whenever needed.<jupyter_code>def function(x):
    return np.cos(x)<jupyter_output><empty_output><jupyter_text>Here I calculate the necessary c and d limits of integration. For Monte Carlo Integration, we want to multiply our ratio by the total area of our region (d-c)(b-a). It is rather easy then to see that d should be the maximum of the function and c should be the minimum. However, what if we test on an all positive interval? Suppose then that d > c > 0. In such case, our integral is all wrong because an integral by definition must be to y = 0. Thankfully, the fix is as easy an setting an 'if' case that makes c = 0 if the minimum of the function is greater than 0. Similary, on an all negative interval, setting d = 0 fixes our integral approximation.<jupyter_code>def get_d(f, a, b, N):
    #create necessary arrays
    x_range = np.linspace(a, b, N)
    y_range = np.zeros(N)
    
    for i in range(N):
        y_range[i] = f(x_range[i])
    d = np.amax(y_range)
    
    if(d < 0):
        d = 0
        return d
    
    return d

def get_c(f, a, b, N):
    x_range = np.linspace(a, b, N)
    y_range = np.zeros(N)
    
    for i in range(N):
        y_range[i] = f(x_range[i])
    c = np.amin(y_range)
    
    if(c > 0):
        c = 0
        return c
    
    return c<jupyter_output><empty_output><jupyter_text>Here I use scipy to create an array of N random float values between boundaries (a,b) and (c,d)<jupyter_code>def x_random_function(a, b, N):
    return random.uniform(a, b, N)
def y_random_function(c, d, N):
    return random.uniform(c, d, N)<jupyter_output><empty_output><jupyter_text>Here these functions sorts the arrays of random x and y into arrays that account for 'over' and 'under' the function. In this case, 'over' the function refers to "misses," and 'under' refers to "hits" by the random number generation; moreover, a hit is between the curve and x-axis, and a miss is anywhere else.Besides that, the functions check along the random number arrays sequentially for anything 'over' or 'under' the function, in which case they will be appended to another array accordingly. The resulst is 4 separate arrays of x and y over and under the function respectively.I chose to break up this sorting into positive and negative portions of the function. I also opted to not use 'else' statements when applicable and opted to be explicity careful with including terms in the various >=, , <, 'if,' and 'elif,' statements. It was extra work, but I'm confident my 'if' and 'elif' statements include all possible numbers. I should mention that I did use 'else' for any error debugging I could've overlooked.I'd also like to mention that this can be done easier with two dimensional arrays. Unfortunately, I didn't know how to use two dimensional arrays until after I wrote this function. Given that this function works, I don't see a need to change it.<jupyter_code>def x_over_function(f, a, b, c, d, N, x_mc_copy, y_mc_copy):
    
    x_mc_over = []
    
    for i in range(N):
        if(f(x_mc_copy[i]) >= 0.0):
            if(y_mc_copy[i] >= f(x_mc_copy[i])):
                x_mc_over.append(x_mc_copy[i])
            #this portion accounts for regions below the x-axis and curve
            elif(y_mc_copy[i] < 0.0):
                x_mc_over.append(x_mc_copy[i])
                
        elif(f(x_mc_copy[i]) < 0.0):
            if(y_mc_copy[i] < f(x_mc_copy[i])):
                x_mc_over.append(x_mc_copy[i])
            #this portion accounts for regions above x-axis and curve
            elif(y_mc_copy[i] >= 0.0):
                x_mc_over.append(x_mc_copy[i])
        else:
            print("something terriblely unprecedented has happened! Please check function x_over")
            return 0
            
    return x_mc_over

def y_over_function(f, a, b, c, d, N, x_mc_copy, y_mc_copy):
    
    y_mc_over = []
    
    for i in range(N):
        if(f(x_mc_copy[i]) >= 0.0):
            if(y_mc_copy[i] >= f(x_mc_copy[i])):
                y_mc_over.append(y_mc_copy[i])
            elif(y_mc_copy[i] < 0.0):
                y_mc_over.append(y_mc_copy[i])
                
        elif(f(x_mc_copy[i]) < 0.0):
            if(y_mc_copy[i] < f(x_mc_copy[i])):
                y_mc_over.append(y_mc_copy[i])
            elif(y_mc_copy[i] >= 0.0):
                y_mc_over.append(y_mc_copy[i])
        else:
            print("something terriblely unprecedented has happened! Please check function y_over")
            return 0
            
    return y_mc_over

def x_under_function(f, a, b, c, d, N, x_mc_copy, y_mc_copy):
    
    x_mc_under = []
    
    for i in range(N):
        if(f(x_mc_copy[i]) >= 0.0):
            if(0.0 <= y_mc_copy[i] < f(x_mc_copy[i])):
                x_mc_under.append(x_mc_copy[i])
                
        elif(f(x_mc_copy[i]) < 0.0):
            if(0.0 > y_mc_copy[i] >= f(x_mc_copy[i])):
                x_mc_under.append(x_mc_copy[i])
            
        else:
            print("something terriblely unprecedented has happened! Please check function x_under")
            return 0
    return x_mc_under

def y_under_function(f, a, b, c, d, N, x_mc_copy, y_mc_copy):
    
    y_mc_under = []
    
    for i in range(N):
        if(f(x_mc_copy[i]) >= 0.0):
            if(0.0 <= y_mc_copy[i] < f(x_mc_copy[i])):
                y_mc_under.append(y_mc_copy[i])
                
        elif(f(x_mc_copy[i]) < 0.0):
            if(0.0 > y_mc_copy[i] >= f(x_mc_copy[i])):
                y_mc_under.append(y_mc_copy[i])
        else:
            print("something terriblely unprecedented has happened! Please check function y_under")
            return 0
        
    return y_mc_under<jupyter_output><empty_output><jupyter_text>I decided that for my error tolerance I will be comparing my Monte Carlo Integration to that of the Trapezoidal Integration; it really only seemed right given we closely studied and practiced this technique earlier in the quarter, so what better comparrison to be made against Monte Carlo?Anyway, Trapezoid_core defines the necessary mathematics of the trapezoid integral. Trapezoid_method creates an x array of equal sizes over the entire interval. Then, a tiny step size is created by using the difference between the second and first terms in this array of x's. From there I sequentially call the Trapezoid_core to add up each individual trapezoid.<jupyter_code>def trapezoid_core(f, x_trap, h):
    return 0.5*h*(f(x_trap + h) + f(x_trap))

def trapezoid_method(f, a, b, N):
    x_trap = np.linspace(a, b, N)
    
    h = x_trap[1] - x_trap[0]
    
    trapezoid_integral = 0.0
    
    for i in range(0, len(x_trap) -1, 1):
        trapezoid_integral += trapezoid_core(f, x_trap[i], h)
        
    return trapezoid_integral<jupyter_output><empty_output><jupyter_text>This is now the bread and butter and by far the hardest part of the program. I initialize an array for the progression of the Monte Carlo Integral over the iterations for a fun plot of error at the end of this program.I first initalize positive, negative, and ambigious totals. This is because I have to be very careful to subtract a Monte Carlo Integral if it's below the curve. Anyway, I create a loop that checks if the random y is greater than 0 and less than the function. In which case, the integral there will be positive so I add a count to the total positive. If the above case isn't satisfied, I check to see if the random y is greater than the function but less than 0. In which case, the integratl there will be negative so I add a count to the total negative. If neither of those are satisfied, I don't add anything to total positive or total negative. It's not necessary to write plus 0, but I think it gets my idea across better then leaving nothing there. Regardless of the above conditions, I must add count to the total.What proceeds that bulk is calculating the positive and negatives integrals and finding their difference; this difference of which is the true Monte Carlo Integral. I append a value to the array of integral values and then check to see if I've hit the tolerance of the trapezoid integral passed into this function. If I've hit the tolerance, there's no need to continue so I can simply return the integral array. If I never hit the tolerance, that's okay maybe I was unlucky or need more iterations.I would like to note that to a brief glance, it appears silly to multiply the total positive by the entire iterval (d-c)(b-a) and likewise for the total negative. However, if one were to do the math, we see that this (d-c)(b-a) can simply be distributed out of the calculation 'mc_integral = mc_integral_positive - mc_integral_negative.' Therefore, the multiplication of (d-c)(b-a), though maybe redundant, is not contributing to error.<jupyter_code>def montecarlo_integration(f, a, b, c, d, N, tolerance, x_mc, y_mc, trapezoid_integral):
    #we will need this for later
    mc_integral_array = []
    
    total_under_positive = 0.0
    total_under_negative = 0.0
    total = 0.0
    
    for i in range(N):
        if(0 <= y_mc[i] < f(x_mc[i])):
            total_under_positive += 1
        elif(0 > y_mc[i] >= f(x_mc[i])):
            total_under_negative += 1
        else:
            total_under_negative += 0
            total_under_positive += 0
            
        total += 1
        
        mc_integral_positive = (total_under_positive/total)*(d-c)*(b-a)
        mc_integral_negative = (total_under_negative/total)*(d-c)*(b-a)
        mc_integral = mc_integral_positive - mc_integral_negative
        mc_integral_array.append(mc_integral)
        
        if(np.abs(mc_integral - trapezoid_integral) <= tolerance):
            print("we've reached our tolerance at", i, "iterations!")
            return mc_integral_array
        
    return mc_integral_array<jupyter_output><empty_output><jupyter_text>Here I declare the boundars of interest (a,b). From there, I use function, get_c, and get_d to calculate the verticals limits of integration. I also declare the number of iterations, typically convering to the tolerance between 20,000 and 80,000 iterations.I will call the random number generators to get our random x and y arrays for the rest of the program. I will also call trapezoid_method to get a trapezoid sum for the rest of the program. Finally, I will call montecarlo_integration to calculate the Monte Carlo Integral and print the last content in the integral array corresponding to the true integral. For comparison I will also print the trapezoid integral.<jupyter_code>#lets set our bounds
N = 100000
a = 0.0
b = 1.75
c = get_c(function, a, b, N)
d = get_d(function, a, b, N)
tolerance = 1.0e-6

randomx = x_random_function(a, b, N)
randomy = y_random_function(c, d, N)

trapezoid_sum = trapezoid_method(function, a, b, N)
mc_sums = montecarlo_integration(function, a, b, c, d, N, tolerance, randomx, randomy, trapezoid_sum)

print("This is the trapezoid method:", trapezoid_sum)
print("This is the Monte Carlo Method:", mc_sums[-1])<jupyter_output><empty_output><jupyter_text>I will now plot the results!I call the x and y over and under functions scatter. I then create a range of x values and call function to plot our function at each point. The x limits are (a,b) and the y limits are (c,d). Scattering, labeling, coloring, legending, and adjusting point sizes accordingly and the plot is finished!I couldn't figure out a pretty way to display the Monte Carlo Integral value over the graph, so under Matt's advice I decided to use figtext which is notoriously tricky. As such, I have these rather bothersome '' and () in the display.<jupyter_code>x_random_misses = x_over_function(function, a, b, c, d, N, randomx, randomy)
y_random_misses = y_over_function(function, a, b, c, d, N, randomx, randomy)
x_random_hits = x_under_function(function, a, b, c, d, N, randomx, randomy)
y_random_hits = y_under_function(function, a, b, c, d, N, randomx, randomy)

x_graph = np.linspace(a, b, N)
y_graph = function(x_graph)
plt.plot(x_graph, function(x_graph), label = 'f(x) = cos(x)', c = 'black')

plt.xlim([a, b])
plt.ylim([c, d])
plt.xlabel('x')
plt.ylabel('f(x)')

plt.scatter(x_random_misses, y_random_misses, c = 'red', s = 0.5, label = 'red points are misses')
plt.scatter(x_random_hits, y_random_hits, c = 'blue', s = 0.5, label = 'blue points are hits')

plt.hlines(0, a, b, color = "black")
plt.legend(loc = 'upper right', framealpha = 1.0)
plt.figtext(0.3, 1, s = f'{"F(a,b) = ", mc_sums[-1]}')<jupyter_output><empty_output><jupyter_text>Cos(x) can actually be solved analytically, so I decided it would be fun to produce a plot of my approximation for the integral over the number of iterations and compare to the analytic solution. Much of the coding here is the same as above, although I should note that I cannot easily plot the trapezoid integral approximation here.The trapezoid_method uses all N, whereas montecarlo_integration uses necessary N up to the tolerance. Therefore, I cannot plot these on the same interval. I could simply recall and redefine the trapezoid_method for the length of the array of integrals, but I don't see a point to that as this trapezoid integral will now be different to the one calculate above.The important takeaway from this fun little plot is that the function montecarlo_integration does appear to converge and my code works! This is proof that I did not get lucky when I called the function as Monte Carlo Integration is a gamblers game.<jupyter_code>#I will now plot i vs mc_sums
analytic_solution = np.sin(7/4)
iterations = np.arange(0, len(mc_sums), 1)

plt.hlines(0, 0, len(mc_sums)-1, color = "black")
plt.hlines(analytic_solution, 0, len(mc_sums)-1, color = "black", label = "analytic solution")


plt.xlim([0, len(mc_sums)-1])
plt.ylim(np.amin(mc_sums), np.amax(mc_sums))
plt.xlabel('iterations i')
plt.ylabel('integral value')

plt.plot(iterations, mc_sums, c = 'red', label = 'mc_sums[i]')

plt.legend(loc = 'lower right')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08623663332183512, 'fraction_numerical': 0.008347706105553639, 'mean_word_length': 3.4235581324382056, 'pattern_counts': {'":': 0, '<': 59, '<?xml version=': 0, '>': 60, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '632226', 'license': 'MIT', 'n_tokens_mistral': 4783, 'n_tokens_neox': 4530, 'n_words': 1796, 'path': 'Final Project/Monte Carlo Integration Final Project  FOR REAL.ipynb', 'repo_name': 'Atlasshrugs00/astr-119', 'size': '14502'}