MIT 6.0002 PS4: Simulating the Spread of Disease and Bacteria Population

请注意,本文编写于 122 天前,最后修改于 122 天前,其中某些信息可能已经过时。

本文是我对 MIT 6.0002 课程(Introduction to Computational Thinking and Data Science)配套作业 Problem Set 4 的解答。这次,我们要完成一个模拟细菌和抗生素相互作用的计算模型。

问题描述

详细的要求见原题,从 MIT OCW 下载原题:Problem Set 4: Simulating the Spread of Disease and Bacteria Population (ZIP) (This ZIP file contains: 1 .pdf file and 2 .py files)

Problem 1: Implementing a Simple Simulation (No Antibiotic Treatment)

在“简单模拟”中,没有抗生素的参与,细菌也不具有抗性。

SimpleBacteria 类

SimpleBacteria的构造函数只需要给成员变量赋值,分别为细菌的出生率和死亡率。

def __init__(self, birth_prob, death_prob):
    self.birth_prob = birth_prob
    self.death_prob = death_prob

is_killed依概率判断细菌是否死亡,返回布尔类型。我们使用random.random()来产生一个 $[0, 1) $ 之间的随机实数。

def is_killed(self):
    return random.random() <= self.death_prob

根据 doc string 的中的说明完成reproduce方法。在每个时间点,依概率判断在当前时刻是否会繁殖新的细菌——如果是,则返回一个新的SimpleBacteria 对象;如果不是,则抛出异常。

def reproduce(self, pop_density):
    if random.random() <= self.birth_prob * (1 - pop_density):
        return SimpleBacteria(self.birth_prob, self.death_prob)
    else:
        raise NoChildException('The cell does not reproduce')

Patient 类

构造函数只需要给成员变量赋值,分别为细菌的列表和最大数量。

def __init__(self, bacteria, max_pop):
    self.bacteria = bacteria
    self.max_pop = max_pop

get_total_pop返回细菌数量。

def get_total_pop(self):
    return len(self.bacteria)

根据update的 doc string 完成该方法。首先创建一个新列表surviving_bacteria,遍历所有细菌,将那些未死亡的放入这个列表;然后根据公式计算细菌密度;使用该密度来判断surviving_bacteria中的每个细菌是否繁殖,将那些新繁殖出来的细菌放入offspring列表;最后,将两个临时列表合并并赋值给self.bacteria。注意调用reproduce的语句要放在try块中,因为其中可能会抛出异常,我们要进行异常处理。

def update(self):
    # 1
    surviving_bacteria = []
    for bacterium in self.bacteria:
        if not bacterium.is_killed():
            surviving_bacteria.append(bacterium)
    
    # 2
    density = len(surviving_bacteria) / self.max_pop

    # 3
    offspring = []
    for bacterium in surviving_bacteria:
        try:
            offspring.append(bacterium.reproduce(density))
        except NoChildException:
            continue

    # 4
    self.bacteria = surviving_bacteria + offspring
    
    return self.get_total_pop()

Problem 2: Running and Analyzing a Simple Simulation (No Antibiotic Treatment)

calc_pop_avg

populations参数是一个矩阵(其实是列表的列表),其每一行代表一次试验(trial),每一列代表一个时刻(time step)。calc_pop_avg函数计算给定时刻的所有试验的平均细菌数量。

def calc_pop_avg(populations, n):
    total_size = 0
    for population in populations:
        total_size += population[n]
    return total_size / len(populations)

或者,本函数还可以使用 Numpy 来实现,如下:

def calc_pop_avg(populations, n):
    return np.array(populations).mean(axis=0)[n]

以上两种实现方式均可通过ps4_tests.py中的单元测试:

test_calc_pop_avg (__main__.ps4_calc) ... 762.5
ok

simulation_without_antibiotic

首先根据参数创建populations空列表(实际上是空矩阵),然后进行num_trials次试验。每次试验中,创建bacteria列表和patient对象用于本次试验,将一次试验的结果保存在population(单数)列表中。一次试验结束后,population加到populations后面,称为矩阵的新一行。

绘制图像的y_coords参数使用列表解析来构建,这里可以使用我们已经实现好了的calc_pop_avg函数来构建这个列表:[calc_pop_avg(populations, step) for step in range(timesteps)]

def simulation_without_antibiotic(num_bacteria,
                                  max_pop,
                                  birth_prob,
                                  death_prob,
                                  num_trials):
    populations = []
    timesteps = 300
    for i in range(num_trials):
        population = []  # 本次trial在每个时间的的细菌数量
        bacteria = [SimpleBacteria(birth_prob, death_prob)] * num_bacteria
        patient = Patient(bacteria, max_pop)

        for step in range(timesteps):
            population.append(patient.update())
            
        populations.append(population)

    # 绘制图像
    make_one_curve_plot(range(timesteps),
        [calc_pop_avg(populations, step) for step in range(timesteps)],
        'Timestep', 'Average Population', 'Without Antibiotic')
    
    return populations

运行模拟

populations = simulation_without_antibiotic(100, 1000, 0.1, 0.025, 50)

模拟完毕后,将运行模拟的语句注释掉。

Problem 3: Calculating a Confidence Interval

根据公式计算期望、标准差和置信区间即可。为了保持代码风格的前后一致,我这里也没有使用 Numpy,而是完全按照公式手动实现了一遍,并不是很复杂。

计算标准差的函数:

def calc_pop_std(populations, t):
    mean = calc_pop_avg(populations, t)  # 期望
    standard_deviation = 0
    for population in populations:
        standard_deviation += (population[t] - mean) ** 2
    standard_deviation = math.sqrt(standard_deviation / len(populations))  # 标准差
    return standard_deviation

95% 置信区间:

def calc_95_ci(populations, t):
    mean = calc_pop_avg(populations, t)
    sem = calc_pop_std(populations, t) / math.sqrt(len(populations))
    return mean, 1.96 * sem

至此,ps4_test.py中的单元测试可以全部通过了:

test_calc_95_ci (__main__.ps4_calc) ... 6.653904117133038
ok
test_calc_pop_avg (__main__.ps4_calc) ... 762.5
ok
test_calc_pop_std (__main__.ps4_calc) ... 10.735455276791944
ok

----------------------------------------------------------------------
Ran 3 tests in 0.005s

OK

我们计算第 299 个时刻的 95% 置信区间:

mean, sem = calc_95_ci(populations, 299)
print('95% confidence interval: {:.3f} ± {:.3f}'.format(mean, sem))

得到结果:

95% confidence interval: 761.520 ± 4.821

Problem 4: Implementing a Simulation with an Antibiotic

这里,考虑抗生素的作用,同时,也要考虑细菌的抗性。

ResistantBacteria 类

不是所有细菌都具有抗性,一个细菌是否具有抗性是在该细菌出生时决定的。细菌有一个参数mut_prob表示突变概率,当细菌突变时,细菌即获得抗性,同时其死亡率也会受到影响。

构造函数:

def __init__(self, birth_prob, death_prob, resistant, mut_prob):
    SimpleBacteria.__init__(self, birth_prob, death_prob)
    self.resistant = resistant
    self.mut_prob = mut_prob

get_resistant

def get_resistant(self):
    return self.resistant

is_killed需要根据细菌是否具有抗性来执行不同的行为,具有抗性的细菌的死亡率更大:

def is_killed(self):
    if self.get_resistant():
        return random.random() <= self.death_prob
    else:
        return random.random() <= self.death_prob / 4

reproduce的 doc string 写得非常详细,没有什么需要补充说明的。

def reproduce(self, pop_density):
    if random.random() <= self.birth_prob * (1 - pop_density):
        if self.get_resistant():
            return ResistantBacteria(self.birth_prob, self.death_prob, self.resistant, self.mut_prob)
        else:
            offspring_resistant = (random.random() <= self.mut_prob * (1-pop_density))  # bool,表示后代是否有抗性
            return ResistantBacteria(self.birth_prob, self.death_prob, offspring_resistant, self.mut_prob)
    else:
        raise NoChildException('The cell does not reproduce')

TreatedPatient 类

该子类有一个新的布尔型成员变量self.on_antibiotic,表示该病人是否使用了抗生素。初始时,该值为 False。

def __init__(self, bacteria, max_pop):
    Patient.__init__(self, bacteria, max_pop)
    self.on_antibiotic = False

set_on_antibiotic方法表示对病人使用抗生素。

    def set_on_antibiotic(self):
        self.on_antibiotic = True

get_resist_pop方法用于获取具有抗性的细菌数量,我使用了列表解析的方法构造出一个由所有带抗性细菌组成的列表,然后返回它的长度。

def get_resist_pop(self):
    return len([bac for bac in self.bacteria if bac.get_resistant()])

根据update的 doc string 完成该方法,下面代码中的注释序号与 doc string 中的序号 1~5 相对应。

def update(self):
    # 1
    surviving_bacteria = []
    for bacterium in self.bacteria:
        if not bacterium.is_killed():
            surviving_bacteria.append(bacterium)
    
    # 2
    if self.on_antibiotic:
        surviving_bacteria = [bac for bac in surviving_bacteria if bac.get_resistant()]

    # 3
    density = len(surviving_bacteria) / self.max_pop

    # 4
    offspring = []
    for bacterium in surviving_bacteria:
        try:
            offspring.append(bacterium.reproduce(density))
        except NoChildException:
            continue

    # 5
    self.bacteria = surviving_bacteria + offspring
    
    return self.get_total_pop()

Problem 5: Running and Analyzing a Simulation with an Antibiotic

由于原题变量命名不规范且混乱,因此以下函数中有些变量可能容易引起误解,现说明:populationsresistant_pop保存的是所有 trial 中所有时刻的信息,它们都是矩阵(列表的列表);population(单数)和res_population是单次试验的信息,它们都是一维列表;每次试验结束后,分别将populationres_population append 到populationsresistant_pop中。

def simulation_with_antibiotic(num_bacteria,
                               max_pop,
                               birth_prob,
                               death_prob,
                               resistant,
                               mut_prob,
                               num_trials):
    populations = []
    resistant_pop = []
    timesteps_1 = 150
    timesteps_2 = 250

    for i in range(num_trials):
        population = []       # 本次trial在每个时间的的细菌数量
        res_population = []   # 本次trial在每个时间的抗性细菌数量

        bacteria = [ResistantBacteria(birth_prob, death_prob, resistant, mut_prob)] * num_bacteria
        patient = TreatedPatient(bacteria, max_pop)

        # 使用抗生素之前
        for step in range(timesteps_1):
            population.append(patient.update())
            res_population.append(patient.get_resist_pop())

        # 使用抗生素之后
        patient.set_on_antibiotic()
        for step in range(timesteps_2):
            population.append(patient.update())
            res_population.append(patient.get_resist_pop())

        # 完成一次trial
        populations.append(population)
        resistant_pop.append(res_population)

    # 绘制图像
    timesteps = timesteps_1 + timesteps_2
    make_two_curve_plot(range(timesteps),
        [calc_pop_avg(populations, step) for step in range(timesteps)],
        [calc_pop_avg(resistant_pop, step) for step in range(timesteps)],
        'Total', 'Resistant', 'Timestep', 'Average Population',
        'With an Antibiotic'
    )

    return populations, resistant_pop

Simulation A

运行第一个模拟,并打印第 299 个时刻的 95% 置信区间:

total_pop, resistant_pop = simulation_with_antibiotic(num_bacteria=100,
                                                      max_pop=1000,
                                                      birth_prob=0.3,
                                                      death_prob=0.2,
                                                      resistant=False,
                                                      mut_prob=0.8,
                                                      num_trials=50)

mean, sem = calc_95_ci(total_pop, 299)
print('95% confidence interval for the total population: {:.3f} ± {:.3f}'.format(mean, sem))
mean, sem = calc_95_ci(resistant_pop, 299)
print('95% confidence interval for the resistant bacteria: {:.3f} ± {:.3f}'.format(mean, sem))

第一个模拟的置信区间:

95% confidence interval for the total population: 202.000 ± 8.705
95% confidence interval for the resistant bacteria: 202.000 ± 8.705

Simulation B

运行第二个模拟,并打印第 299 个时刻的 95% 置信区间:

total_pop, resistant_pop = simulation_with_antibiotic(num_bacteria=100,
                                                      max_pop=1000,
                                                      birth_prob=0.17,
                                                      death_prob=0.2,
                                                      resistant=False,
                                                      mut_prob=0.8,
                                                      num_trials=50)

mean, sem = calc_95_ci(total_pop, 299)
print('95% confidence interval for the total population: {:.3f} ± {:.3f}'.format(mean, sem))
mean, sem = calc_95_ci(resistant_pop, 299)
print('95% confidence interval for the resistant bacteria: {:.3f} ± {:.3f}'.format(mean, sem))

第二个模拟的置信区间:

95% confidence interval for the total population: 0.000 ± 0.000
95% confidence interval for the resistant bacteria: 0.000 ± 0.000

Problem 6 Write-up

其实上面已经放过图像了,这里按照题目要求再放一次。下面三幅图和上文中的图像完全相同。

Problem 2 的图像

Problem 5 的图像

  1. 在 Simulation A 中,对患者使用抗生素之前,细菌总数量逐渐增加,并在第 100 个时刻后接近max_pop。在 Simulation B 中,细菌总数量先增加后减少,这是因为 B 的出生率比 A 要低近一半。
  2. 在 Simulation A 和 B 中,抗性细菌数量的变化趋势都是先增加后减少,只是增减的幅度略有不同,这是由于 A 和 B 中细菌具有不同的出生率,A 的出生率高,因此抗性细菌数量增长得更快,而且达到的最大值也比 B 中的更大。
  3. 使用抗生素后,细菌数量骤减,非抗性细菌全部死亡,因此使用抗生素后 total 等于 resistant。在 Simulation A 中,出生率足够高,因此抗性细菌在抗生素的作用下仍能增加,且渐进于一个稳定的值;在 Simulation B 中,出生率较低,所以抗性细菌几乎全部死亡,难以繁殖。
  4. 使用抗生素后,仍能存活的细菌只有抗性细菌。根据细菌的出生率的不同,细菌种群数量的变化趋势也有所不同,这一点在 3 中已经说明。

总结

  • 计算模型的概念。本次实验使用 Python 完成了一个模拟细菌和抗生素相互作用的计算模型,计算模型(computational model)是指使用大量的计算资源来用计算机模拟研究一个复杂系统的行为。
  • 伪随机数和种子。熟悉任何一种编程语言的人都知道伪随机数的概念——伪随机数是用确定性的算法计算出来的均匀分布的随机数序列,并不真正的随机,但具有类似于随机数的统计特征。一个随之关联的概念称为“种子(seed)”,对于一个相同的 seed,多次运行同一段程序,得到的随机序列是相同的。这个特性有利于我们 调试程序。在 Python 中设置伪随机数种子的语句是random.seed(0)。调试完成后,应该删除这条语句。
  • 异常处理。Python 的异常处理使用的是try/except语句,该语句检测try中的错误(即抛出的异常),然后用except语句捕获异常并给予相应的处理。在本次实验中,except语句中仅有一个continue。在生产环境中,情况将会复杂得多。
  • 函数局部变量与参数重名导致难以发现的错误。Python 中使用一个变量前并不需要声明它(至少不需要像 C/C++ 那样显式地声明),我一直以为这是一个优点,直到今天,我在一个函数中错误地使用了一个变量名,而这个变量名恰好是该函数的一个参数,但我并没有意识到这一点,以为它仅是一个普通的局部变量而已。调试了很久才发现这个错误。

Comments

添加新评论

已有 4 条评论

好像也没什么问题。。。

别吧,评论还要审核?
# 注释后的代码颜色有些问题

Jed Jed 回复 @爸爸

为了防止垃圾评论。。。一个邮箱审核一次就可以了

沙发,学习了!赞! ୧(๑•̀⌄•́๑)૭