Python蒙特卡洛模拟结果绘图异常：优化代码速度快但图形错误

原代码（运行正确但耗时久） #

以下代码可正确绘制Monte Carlo模拟结果，但运行耗时较长，生成的图形符合预期：

runs = 1000000

import matplotlib.pyplot as plt
import random
import numpy as np
from scipy.stats import norm
from scipy.stats import uniform
import seaborn as sns
import pandas 

def steadystate():
    p=0.88
    Cout=4700000000
    LambdaAER=0.72
    Vol=44.5
    Depo=0.42
    Uptime=0.1
    Effic=0.38
    Recirc=4.3
    x = random.randint(86900000,2230000000000)
    conc = ((p*Cout*LambdaAER)+(x/Vol))/(LambdaAER+Depo+(Uptime*Effic*Recirc))
    return conc

x = 0
while x < runs:
    #results = steadystate (Faster)
    results = np.array([steadystate() for _ in range(runs)])
    print(results)
    x+=1
    
ax = sns.distplot(results,
                  bins=100,
                  kde=True,
                  color='skyblue',
                  hist_kws={"linewidth": 15,'alpha':1})
ax.set(xlabel='Uniform Distribution ', ylabel='Frequency')

优化后代码（速度提升但绘图异常） #

以下代码运行速度大幅提升，但无法正确生成分布图，仅输出异常图形：

import numpy.random
runs = 1000000000

import matplotlib.pyplot as plt
import random
import numpy as np
from scipy.stats import norm
from scipy.stats import uniform
import seaborn as sns
import pandas 

#Simulation Function
def steadystate(count):  
    p=0.88
    Cout=4700000000
    LambdaAER=0.72
    Vol=44.5
    Depo=0.42
    Uptime=0.1
    Effic=0.38
    Recirc=4.3
    x = numpy.random.randint(86900000, 2230000000000, dtype=np.int64)  

    conc = ((p*Cout*LambdaAER)+(x/Vol))/(LambdaAER+Depo+(Uptime*Effic*Recirc))
    return conc

x = 0
while x < runs:
    results = steadystate(runs)
    print(results)
    x+=1
#While loop to run iterations of the simulation for x number of runs
    
#Using seaborn to plot the results in a uniform distribution
ax = sns.distplot(results,
                  bins=100,
                  kde=True,
                  color='skyblue',
                  hist_kws={"linewidth": 15,'alpha':1})
ax.set(xlabel='Uniform Distribution ', ylabel='Frequency')

问题诊断与解决建议 #

核心问题 #

1. steadystate函数逻辑错误：函数接收count参数但未使用，numpy.random.randint未指定size参数，每次仅生成单个随机数，返回的是标量而非数组。
2. 外层while循环完全冗余且错误：runs设为1e9，循环会执行10亿次，每次覆盖results变量，最终results仅保留最后一次生成的单个数值，导致Seaborn无法绘制分布直方图，只能显示单个点的异常图形。
3. 模块导入冗余混乱：重复导入numpy相关模块，无意义。

修正后的代码 #

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

runs = 1000000  # 建议根据内存调整，1e9会占用极大内存，原代码用1e6更合理

def steadystate(count):  
    p = 0.88
    Cout = 4700000000
    LambdaAER = 0.72
    Vol = 44.5
    Depo = 0.42
    Uptime = 0.1
    Effic = 0.38
    Recirc = 4.3
    # 生成count个随机数组成的数组
    x = np.random.randint(86900000, 2230000000000, size=count, dtype=np.int64)  
    # 向量运算，直接生成对应数量的结果数组
    conc = ((p * Cout * LambdaAER) + (x / Vol)) / (LambdaAER + Depo + (Uptime * Effic * Recirc))
    return conc

# 直接生成所有模拟结果，无需循环
results = steadystate(runs)

# 绘制分布图
ax = sns.distplot(results,
                  bins=100,
                  kde=True,
                  color='skyblue',
                  hist_kws={"linewidth": 1, 'alpha': 0.8})  # 调整linewidth避免图形过粗
ax.set(xlabel='Uniform Distribution', ylabel='Frequency')
plt.show()  # 需显式调用show()显示图形

关键优化点 #

• 利用numpy的向量运算一次性生成所有随机数和结果，避免Python循环的开销，同时保证返回的是数组类型。
• 移除无意义的while循环，直接调用一次函数即可获取所有模拟结果。
• 调整hist_kws中的linewidth参数，避免原代码中15的线宽导致图形显示异常。
• 添加plt.show()确保图形正常显示。
• 注意：runs=1e9会占用极大内存（单个int64是8字节，1e9就是8GB），建议根据实际内存情况调整为合理数值（如原代码的1e6）。

内容的提问来源于stack exchange，提问作者bowlcutinmatrix