在Python中模拟的衰变链

Question

我正在尝试制作一个计算机程序,可以模拟氡-222的1,000,000个原子的衰变(进入Pol -ium-218进入Lead-214进入铋-214进入铅-210).每个原子在一个单位时间后有0.000126的机会腐烂成Pol -218.我有一个1,000,000 0 m的列表分配给变量Rn来表示Radon-222.然后,为了确定特定元素是否衰减,我选择1到1,000,000之间的随机整数:如果此数字小于126,则Rn [j]的值从0切换为1(衰减),并且我们将0添加到在单位时间之后,我预计会有大约126个Pol原子,在另一个单位时间之后,我们应该有大约126个额外的Pol原子,但是其中一些可能会腐烂到Lead-214中. .

我们将时间视为此建模中的离散实体.在此处查看10个单位时间后的理想值列表.

这是一个有效的代码,但速度很慢.

import random

def one(s):
    count=0
    for i in range(len(s)):
        if s[i]==1:
            count+=1
    return count

Rn=[1]*1000000
Po=[0]*1000000
Pb214=[0]*1000000
Bi=[0]*1000000
Pb210=[0]*1000000

print([0,1000000,0,0,0,0])

for i in range(1,100):
    for j in range(len(Rn)):
        rndecay=random.random()
        if rndecay<=0.000126:
            Rn[j]=2        

    for j in range(len(Po)):
        if Po[j]==1:
            podecay=random.random()
            if podecay<=0.203:
                Po[j]=2

        if Rn[j]==2 and Po[j]==0:
            Po[j]=1

    for j in range(len(Pb214)):            
        if Pb214[j]==1:
            decay214=random.random()
            if decay214<=0.0255:
                Pb214[j]=2
        if Po[j]==2 and Pb214[j]==0:
            Pb214[j]=1

    for j in range(len(Bi)):
        if Bi[j]==1:
            bidecay=random.random()
            if bidecay<=0.0346:
                Bi[j]=2
        if Pb214[j]==2 and Bi[j]==0:
            Bi[j]=1

    for j in range(len(Pb210)):
        if Bi[j]==2 and Pb210[j]==0:
            Pb210[j]=1

    print([i,one(Rn),one(Po),one(Pb214),one(Bi),one(Pb210)])

我该如何优化此代码？

Answer 1

我认为你可能会以不同的方式更好地解决问题.您可以只跟踪每个物种的原子数,而不是拥有1,000,000个元素列表吗？

具有衰减概率p的n群体中的衰变数遵循具有这些参数的二项分布.许多基本的统计软件包将允许您生成遵循分布的随机数(而不是使用的统一分布random.random).因此,您只需为每个可能的转换提取一个随机值,然后使用结果更新计数.这种方法完全等同于原始帖子中使用的方法.

这是一个使用我编制的具有转移概率的计数的示例:

from numpy.random import binomial

atoms = {'Rn222': 1000000,
         'Po218': 0,
         'Pb214': 0,
         'Bi214': 0,
         'Pb210': 0}

for _ in range(100):
    Rn222_Po218 = binomial(atoms['Rn222'], 0.000126)
    Po218_Pb214 = binomial(atoms['Po218'], 0.001240)
    Pb214_Bi214 = binomial(atoms['Pb214'], 0.003450)
    Bi214_Pb210 = binomial(atoms['Bi214'], 0.012046)

    atoms['Rn222'] -= Rn222_Po218
    atoms['Po218'] += Rn222_Po218

    atoms['Po218'] -= Po218_Pb214
    atoms['Pb214'] += Po218_Pb214

    atoms['Pb214'] -= Pb214_Bi214
    atoms['Bi214'] += Pb214_Bi214

    atoms['Bi214'] -= Bi214_Pb210
    atoms['Pb210'] += Bi214_Pb210

print atoms

编辑添加示例.编辑从评论中添加二项式解释.