Build n distribution of Monte Carlo Box Method in Python

user_9952960 0
Viewed 103

I'm trying to perform resolve this exercise:

enter image description here

using the Monte Carlo Method "Box Method".

The code that i implemented is this:

#UNIDIMENSIONAL INTEGRATION

import numpy as np                            #library for numerical calculations
import matplotlib.pyplot as plt               #library for plotting purposes
from scipy import random                      #needed for generate random number
from sympy import symbols, integrate, exp     #needed for integrate function
from scipy.stats import norm                  #needed for gaussian fit

#*******************************************************************************

def f(x,n):                              #definition of the function to integrate

  return x**n

#*******************************************************************************

for i in range(1, 6):                    #for cycle over the period

  x = symbols('x')                       #needed for the integration

  print("The exact mathematical value of the integral with dimension N", i, "is:", integrate(f(x,i),(x, 0,1)).evalf(2), "\n") 

print("************************************************************************* \n")

#*******************************************************************************

N = 10**3                                 #number of point generated, statistics


for j in range(1,6):                      #for loop over the D dimensions

  ans = 0                                 #variable ans 

  list_ans = []                           #list to store all the values for plotting. 

  n_below_curve = 0                       #variable n_below_curve

  for k in range(N):

    for i in range(N):

      x0 = np.random.uniform(0,1)
      y0 = np.random.uniform(0,j)  

      if (y0 <= f(x0,j)):

        n_below_curve += 1

    ans = (n_below_curve/N) * (1*j)

    list_ans.append(ans)    
    
  print("\nThe distribution of the results of integral with dimension N", j, ".\n")

  _, bins, _ = plt.hist(plt_vals, int(np.sqrt(N)), density=True)          #sintex to create a histogram from a dataset x with n bins 
                                                                          #and store an array specifying the bin ranges in the variable bins.

  mu, sigma = norm.fit(plt_vals)                                          #get the mean and standard deviation of data

  best_fit_line = norm.pdf(bins, mu, sigma)                               #get a line of best fit for the data

  print("\n")

  print("The distribution of the results of integral with dimension N", n, ".\n")

  print("The mean of the distribution is ", mu, ". The sigma of the distribution is", sigma, ".\n") 

  plt.hist(plt_vals, bins, ec="black", density=True)                       #compute and draw the histogram of x with n bins
  plt.plot(bins, best_fit_line)                                            #plot y versus x as lines and/or markers
  plt.grid()                                                               #configure the grid lines
  plt.xlabel("Results of integral")                                        #set the label for the x-axis 
  plt.ylabel("N")                                                          #set the label for the y-axis 
  plt.title("Distribution of results of integral")                         #set a title for the histogram
  plt.show()                                                               #display all open figures

#*******************************************************************************

But the outputs give 5 distribution with only one bin not empty...

enter image description here

Why this happened?

The code without the distribution part works perfectly, so I suppose that the problem arise when I try to build the histograms...

Thanks in advance.

montecarlo python
1 Answers

Not sure why do you want to get points below the curve. You already have integral from 0 to 1, just compute the mean value of the function

Code, Python 3.9 Win 10x64

import numpy as np

rng = np.random.default_rng()

N = 100000

acc = np.zeros(6) # accumulator

def fn(x, n):
    return np.power(x, n)

for n in range(1, 6):
    acc[n] = np.mean(fn(rng.random(N), n))

for n in range(1, 6):
    print((n, acc[n], 1.0/(1+n)))

produced

(1, 0.5020245718854783, 0.5)
(2, 0.3342433303484079, 0.3333333333333333)
(3, 0.24902862198284664, 0.25)
(4, 0.20002409001063104, 0.2)
(5, 0.16627446402760915, 0.16666666666666666)
user_40446960
Related