python - 使用 NumPy 和 matplotlib 绘制的图形上未显示某些点

Question

对于以下代码，其工作是对函数 f 执行蒙特卡洛积分，我想知道如果将 f 定义为 y = sqrt(1-x^2)(这是单位四分之一的方程)会发生什么圆，并指定一个大于 1 的端点，因为我们知道 f 仅为 0

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)

令我惊讶的是，该程序仍然有效，并给了我下面的图片。

我怀疑它有效，因为在 NumPy 中，在对数组执行操作时，数组中的 nan 会被忽略。但是，我不明白为什么图片只包含 x 和 y 坐标都在 0 到 1 之间的点。不在这个范围内但其值通过以下方式计算的点在哪里

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

Answer 1

您可以打印出数组(例如，仅生成一个随机点)，然后看到它们不会进入 ind_below也不ind_above ...

这是因为所有涉及 nan 的比较返回False 。 (另请参阅:What is the rationale for all comparisons returning false for IEEE754 NaN values?)。 (因此 y_rand < nan 和 y_rand >= nan 的计算结果均为 False )

更改代码的最简单方法是

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

(可选地只计算一次数组)