使用Python计算单次蒙特卡罗模拟中凝聚事件的归一化方差计算

Question

为了在python中使用montecarlo技术，我尝试复制下面的论文（https://pennstate.pure.elsevier.com/en/publications/constant-number-monte-carlo-simulation-of-population-balances）的结果

我有一个问题，了解方差计算应该是什么，或者如果我是完全错误的公式，我正在使用

实际上，问题是计算凝结事件后系统的平均质量（两个粒子合并成一个新粒子，总质量为合并粒子的总和），以蒙特卡罗方式定义接受规则（生成随机数，仅当随机数小于以某种方式定义的概率时才接受凝聚）

然后我要计算标准化方差，绘图是时间的函数，这意味着在模拟过程中要计算几次。我的函数如下所示：

def montecarlorun(N,totalruns,model):

    particles={}
    M_k_hat=[] #number average mass or sample mean
    M_k_bar=[] #mean size
    time=[]
    distributionM=[]

    #Particle initialization of mass
    for i in range(N):
        particles[i]=mass

    #Matrices initialization

    M_k_hat.append(mass) #initial mean mass of the system
    M_k_bar.append(mass)
    time.append(0)


    #Monte Carlo runs

    #variable definitions
    c=0
    k=0
    sigma=[]
    sigmatimes=[]
    pairs=[]
    part=[]
    mean=[]
    meansq=[]
    muolist=[]
    mu2list=[]
    particlesmeanmass=0
    particlesmeanmasssq=0
    kmax=1

    #simulation
    while c<totalruns:
        prtcl1=random.randrange(N)
        prtcl2=random.randrange(N)
        r=random.uniform(0,1) #random number for acceptance probability condition
        if model==1:
            p=1
        else: # definition of probability for a second model, irrelevant to the post as it is
            p,kmax=brownian_kernel(particles[prtcl1], particles[prtcl2],kmax)
        if(r<p):
            newmass=particles[prtcl1]+particles[prtcl2] #newmass after coagulation
            prtcl_copy=random.randrange(N) #copy a particle from the list
            particles.update({prtcl1:newmass}) #update particle 1 with the new mass after the event of coagulation
            copymass=particles[prtcl_copy] #take mass of copied
            particles.update({prtcl2:copymass}) #replace particle 2 with the copy
            time.append((2/(kmax*N))*(N/(N-1))**k) #get time for event
            k=k+1 # event counter
            M_k_hat.append(M_k_hat[k-1]+(copymass/N)) # new mass of system
            M_k_bar.append((N/(N-1))**k) #same
        print(len(M_k_hat))
        print("c=",c)
        #Sigma part new:
        if(c%2000)==0:
            meanmassnow=M_k_hat[-1]
            sigmatimes.append(time[-1])
            for i in range(len(particles)):
                particlesmeanmass=particlesmeanmass+particles[i]/meanmassnow
                particlesmeanmasssq=particlesmeanmasssq+((particles[i]/meanmassnow)**2)
            muo=particlesmeanmass/N
            mu2=particlesmeanmasssq/N
            muolist.append(muo)
            mu2list.append(mu2)
        #Histograms for distributions of masses:
        if(c%2000==0):
            meanmassnow=M_k_hat[-1]
            masses=particles.values()
            masses=[x/meanmassnow for x in masses]
            d=np.histogram(masses,10,density=True)
            d=(d[0] ,0.5*(d[1][1:]+d[1][:-1]))
            distributionM.append([d,time[-1]])

        c=c+1


    for i in range(len(muolist)):
        sigma.append( (mu2list[i]/muolist[i]) - 1 )


    return particles,M_k_bar,M_k_hat,time,sigma,sigmatimes,pairs,distributionM

根据我对论文的理解，我试图从分布的矩计算归一化方差

从论文中，也有一个分析关系。本文的分析结果与模拟结果基本一致。相反，我得到：

其中红色是理论线，蓝色是我的模拟结果。更好地理解蓝色曲线的行为如下：

所以，问题是是否有一个算法问题，或是否有一个错误的数量，我正在绘制（我应该绘制其他东西）。注意，在我看来，从分布的时刻开始，第一时刻总是1