Как правильно обобщить алгоритм Gillesp ie (в идеале объектно-ориентированным способом)? - PullRequest
Новая
       100

Как правильно обобщить алгоритм Gillesp ie (в идеале объектно-ориентированным способом)?

0 голосов
/ 21 января 2020

В настоящее время я работаю над реализацией алгоритма Gillesp ie. Для этого я сначала попытался обобщить код со страницы википедии на модели SIR для подтверждения концепции. Я попытался построить код так, чтобы потом можно было обобщить его таким образом, чтобы не у каждого человека были одинаковые шансы на выздоровление / заражение. Я определил простой класс Person и создал группу из N человек, где каждый человек может быть заражен, восприимчив или вылечен. Затем в основном l oop я создаю кумулятивную сумму вероятностей для «реакции» для каждого человека. После этого я случайным образом выбираю человека для заражения / восстановления и рассчитываю следующий шаг времени. 

import numpy as np

##### CONSTANTS
N = 350         # int; total population
T = 100.0       # float; maximum elapsed time
t = 0.0         # float; start time
_alpha = 10.0   # float; rate of infection after contact
_beta = 0.5     # float; rate of cure
n_I = 1         # int; initial infected population

class Person(object):
    def __init__(self, status="susceptible"):
        super(Person, self).__init__()
        self.status = status
        self.rate = 0
        self.setStatus(status)

    def setStatus(self, status, n=1):
        self.status = status
        if status == "infected":
            self.rate = _beta  # float; rate of cure
        elif status == "susceptible":
            self.rate = _alpha # float; rate of infection after contact
        elif status == "recovered":
            self.rate = 0

persons = [Person() for i in range(N)]
persons[0].setStatus("infected")

t = 0.0         # float; start time
n_I = 1         # int; initial infected population

# Compute susceptible population, set recovered to zero
n_S = N - n_I
n_R = 0

# Main loop
while t < T:
    if n_I == 0:
        break
    R = np.cumsum([person.rate for person in persons])
    u1 = np.random.random()
    R_N = R[-1]

    for i in range(len(R)):
        if R[i] > u1*R_N:
            selected_idx.append(i)
            if persons[i].status == "infected":
                persons[i].setStatus("recovered")
                n_I -= 1
                n_R += 1
            else:
                persons[i].setStatus("infected")
                n_I += 1
                n_S -= 1
            break

    u2 = np.random.random()
    dt = -np.log(u2)/R_N
    t += dt

Однако проблема, с которой я сталкиваюсь, заключается в том, что приведенный выше код не дает таких же результатов, которые я получаю из кода на странице Википедии (см. Результаты моделирования ). Не могли бы вы указать, где я мог бы неправильно понять алгоритм или как я мог изменить свою версию, чтобы получить те же результаты, что и на странице Википедии?
(Полагаю, моя главная проблема - это понимание того, как выбрать правильные скорости реакции?)

Спасибо всем большое.

...