python - 用python计算凝血事件的蒙特卡罗实现中的归一化方差

Question

我正在尝试复制以下论文 ( https://pennstate.pure.elsevier.com/en/publications/constant-number-monte-carlo-simulation-of-population-balances ) 的结果以使用 MONTE Carlo 技术在蟒蛇。

我很难理解方差计算应该是什么，或者我对我使用的方程完全错误。

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

然后我必须计算归一化方差，并且绘图是时间的函数，这意味着在模拟过程中要计算几次。我的功能如下所示：

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

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

从论文中，也有一个解析关系。在本文中，分析结果和模拟结果几乎相同。相反，我得到：

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

所以，问题是是否存在算法问题，或者我绘制的数量是否有错误（我应该绘制其他东西）。请注意，正如我所看到的，从分发的那一刻起，第一个时刻总是一个。

标签： pythonstatisticssimulationmontecarlo