Some points are not displayed on the graph plotted using NumPy and matplotlib

Question

For the following code whose job is to perform Monte Carlo integration for a function f, I was wondering what would happen if I define f as y = sqrt(1-x^2), which is the equation for a unit quarter circle, and specify an endpoint that is greater than 1, since we know that f is only defined for 0<x<1.

import numpy as np
import matplotlib.pyplot as plt

def definite_integral_show(f, x0, x1, N):
    """Approximate the definite integral of f(x)dx between x0 and x1 using
    N random points
    
    Arguments:
    f -- a function of one real variable, must be nonnegative on [x0, x1]
    N -- the number of random points to use
    
    
    """
    #First, let's compute fmax. We do that by evaluating f(x) on a grid
    #of points between x0 and x1
    #This assumes that f is generally smooth. If it's not, we're in trouble!
    x = np.arange(x0, x1, 0.01)
    
    y = f(x)
    print(y)
    f_max = max(y)
    
    
    #Now, let's generate the random points. The x's should be between
    #x0 and x1, so we first create points beterrm 0 and (x1-x0), and 
    #then add x0
    #The y's should be between 0 and fmax
    #
    #                  0...(x1-x0)
    x_rand = x0 + np.random.random(N)*(x1-x0)
    print(x_rand)
    
    y_rand = 0 +  np.random.random(N)*f_max
    
    
    
    #Now, let's find the indices of the poitns above and below
    #the curve. That is, for points below the curve, let's find
    #   i s.t. y_rand[i] < f(x_rand)[i]
    #And for points above the curve, find
    #   i s.t. y_rand[i] >= f(x_rand)[i]
    ind_below = np.where(y_rand < f(x_rand))
    ind_above = np.where(y_rand >= f(x_rand))
    
    
    #Finally, let's display the results
    plt.plot(x, y, color = "red")
    pts_below = plt.scatter(x_rand[ind_below[0]], y_rand[ind_below[0]], color = "green")
    pts_above = plt.scatter(x_rand[ind_above[0]], y_rand[ind_above[0]], color = "blue")
    plt.legend((pts_below, pts_above),
            ('Pts below the curve', 'Pts above the curve'),
            loc='lower left',
            ncol=3,
            fontsize=8)

def f1(x):
    return np.sqrt(1-x**2)

definite_integral_show(f1, 0, 6, 200)

To my surprise, the program still works and gives me the following picture.

I suspect that it works because in NumPy, nan's in an array are just ignored when performing operations on the array. However, I don't understand why the picture only contains points whose x and y coordinates are both between 0 to 1. Where are the points that aren't within this range, but whose values are computed by

x_rand = x0 + np.random.random(N)*(x1-x0)
y_rand = 0 +  np.random.random(N)*f_max

score 2 · Accepted Answer · answered Dec 31 '20 at 03:18

You can just print out the arrays (for example by generating only one random point) and see that they go into neither ind_below nor ind_above...

That's because all comparisons that involves nan returns False. (See also: What is the rationale for all comparisons returning false for IEEE754 NaN values?). (so y_rand < nan and y_rand >= nan both evaluates to False)

The easiest way to change the code is

ind_below = np.where(y_rand < f(x_rand))
ind_above = np.where(~(y_rand < f(x_rand)))

(optionally only compute the array once)

Some points are not displayed on the graph plotted using NumPy and matplotlib

1 Answers1