Python实现蒙特卡洛盒子法积分结果分布时仅单个箱非空的问题排查

Let's break down why your histogram is only showing one non-empty bin—there are a few key issues in your code's loop structure and variable handling that are causing this behavior:

1. Undefined plt_vals Variable #

You’re referencing plt_vals in your histogram code, but this variable is never defined in your script. You clearly meant to use list_ans (the list storing your integral estimates). Every instance of plt_vals needs to be replaced with list_ans.

2. Broken Loop Structure & Unreset Counter #

Your nested loops are making integral estimates accumulate instead of being independent trials:

• You don’t reset n_below_curve at the start of each outer k loop. This means the counter keeps growing with every iteration, turning each ans value into a cumulative estimate (using more and more points) instead of a standalone Monte Carlo trial.
• The result? list_ans fills with values that converge tightly toward the true integral, so almost all estimates cluster around the exact value—hence the single crowded histogram bin.

You need to restructure loops to run independent trials: reset the counter for each trial, generate N points, compute one estimate, and add it to your list.

3. Inefficient Y-Range for the Box Method #

For f(x,j) = x^j on x ∈ [0,1], the maximum function value is 1 (since 1^j = 1). Setting y0 = np.random.uniform(0,j) creates a box taller than necessary (especially for j > 1), which isn’t mathematically wrong but reduces the number of points under the curve and increases variance. Using y0 = np.random.uniform(0,1) is far more efficient.

Corrected Code #

Here’s the fixed version of your script with all issues addressed:

# UNIDIMENSIONAL INTEGRATION
import numpy as np
import matplotlib.pyplot as plt
from scipy import random
from sympy import symbols, integrate, exp
from scipy.stats import norm

def f(x,n):
    return x**n

# Print exact integral values
for i in range(1, 6):
    x = symbols('x')
    exact_val = integrate(f(x,i), (x, 0,1)).evalf(6)
    print(f"The exact mathematical value of the integral for exponent N={i} is: {exact_val}\n")
    print("*************************************************************************\n")

N_POINTS_PER_TRIAL = 10**3  # Number of random points per Monte Carlo trial
N_TRIALS = 10**3            # Number of independent trials to build the distribution

for j in range(1,6):
    list_ans = []
    for _ in range(N_TRIALS):
        n_below_curve = 0
        # Generate points for one trial
        for __ in range(N_POINTS_PER_TRIAL):
            x0 = np.random.uniform(0, 1)
            y0 = np.random.uniform(0, 1)  # Match function's max value on [0,1]
            if y0 <= f(x0, j):
                n_below_curve += 1
        # Calculate integral estimate for this trial
        area_box = 1 * 1  # Area of the bounding box (x: 0-1, y:0-1)
        ans = (n_below_curve / N_POINTS_PER_TRIAL) * area_box
        list_ans.append(ans)
    
    # Plot histogram and Gaussian fit
    print(f"\nDistribution of integral results for exponent N={j}:\n")
    _, bins, _ = plt.hist(list_ans, int(np.sqrt(N_TRIALS)), density=True, ec="black")
    mu, sigma = norm.fit(list_ans)
    best_fit_line = norm.pdf(bins, mu, sigma)
    
    print(f"Mean of distribution: {mu:.6f}, Standard deviation: {sigma:.6f}\n")
    
    plt.plot(bins, best_fit_line, 'r-')
    plt.grid()
    plt.xlabel("Integral Estimate")
    plt.ylabel("Probability Density")
    plt.title(f"Distribution of Monte Carlo Integral Results (N={j})")
    plt.show()

Key Fixes Summary: #

• Replaced undefined plt_vals with list_ans
• Restructured loops to run independent trials (resetting n_below_curve each time)
• Set y-range to [0,1] to match the function's maximum value on the interval
• Renamed variables for clarity (separated points per trial from number of trials)
• Fixed the repeated print statement that referenced an undefined n variable

With these changes, your histograms will show the expected normal distribution of Monte Carlo estimates—centered around the exact integral value, with a spread that narrows as you increase the number of points per trial.

内容的提问来源于stack exchange，提问作者J.Snowden