连续时间SIS模型与算法
易感-感染-易感模型，这个模型最早在 1932 年由 Kermack 和 McKendrick 提出。一般用来描述一些治愈后可再次感染的疾病，如流感与肺结核等疾病。 SIS模型的传播过程可以表述为：易感者(S)通过接触其感染者(I)邻居感染疾病，而感染疾病后的患者会成为新的感染源，同时感染者也会随着时间的推移治愈疾病成为一个易感者。其连续时间模型可以见下图：
其中β为感染速率，μ为恢复速率。对于连续时间的SIS模型，通常采用Gillespie算法，即遗传算法来进行模拟，该算法是基于独立发生的事件，这意味着动力学的每次演变都是由事件决定的。而事件的时间连续性由系统状态的保持时间 𝜏 所分割。具体来说：
(1) 在任意时刻 𝑡，我们根据计算每个个体的转移率 𝜆𝑖(𝑡)，即状态为 S 和状态为 𝐼 的个体， 𝜆𝑖(𝑡) 分别等于 βm 和 μ ，m 表示 i 节点周围的感染者邻居数量。
(2) 然后计算出系统的总速率 𝜔(𝑡) = ∑𝑖 𝜆𝑖(𝑡)，其中系统状态保持时间 𝜏 = 1/𝜔(𝑡)。
(3) 在下一时刻 𝑡 + 𝜏 时，以概率 𝜆𝑖(𝑡)/𝜔(𝑡) 采样选择节点发生状态改变。
(4) 重复 (1) 到 (3) 过程，直到整个系统达到稳定状态或系统中 𝐼 个体为零时停止迭代。
C++代码实现
注：下面代码中网络是直接在模拟中生成的，采用的ER随机网络。通常不建议这样使用，因为每次进行模拟，系统的网络结构就发生了变化。一般来说都是将网络单独生成，再通过导入的方法放入动力学程序当中。
#include <iostream>
#include <ctime>
#include <stdlib.h>
#include <fstream>
#include <cmath>
#include "mt19937ar.c" //从我文章中ER网络实现获取
#include <algorithm>
#include <vector>
using namespace std;
int main()
{
int Network_Size = 1000;
double r = 0.4;
double g = 0.6;
int I0 = 100;
double tend = 100;
double counter = 10;
vector<int> degree(Network_Size,0);
vector<vector<int> > adjlist(Network_Size);
unsigned long idum;
idum=(unsigned long)time(NULL);
init_genrand(idum);
ofstream outfile_ER_ext;
outfile_ER_ext.open("ER_GA_ext.txt");
double average_degree = 4;
int edge = 0.5*Network_Size*average_degree;
for(;;) //生成ER随机网络
{
double judge = 0;
int random_1 = genrand_int32()%Network_Size;
int random_2 = genrand_int32()%Network_Size;
if(random_1 != random_2)
{
for(int i = 0; i < adjlist[random_1].size(); i++)
{
if(adjlist[random_1][i] == random_2)
{
judge = 1;
break;
}
}
if(judge == 1)
{
continue;
}
else
{
adjlist[random_1].push_back(random_2);
degree[random_1]++;
adjlist[random_2].push_back(random_1);
degree[random_2]++;
edge--;
}
}
if(edge == 0)
break;
}
vector<vector<int> > X_vec; //记录所有独立模拟的时间序列
int n_intervals = tend;
double interval = tend/(1.0*n_intervals); //划分取值区间
vector<int> X_vec_ind(n_intervals + 1); //记录每一次的时间序列
vector<double> rate(Network_Size,0); //节点速率
vector<int> I_node(Network_Size,0); //节点状态
vector<double> sum(Network_Size,0);
vector<double> nk(Network_Size,0); //归一化概率
vector<int> kinf(Network_Size,0); //每个节点的感染者邻居数量
for(int n = 0; n < counter; n++) //进行多次独立模拟，用来取平均值
{
double wt = 0;
double dt;
double t = 0;
double I;
int random_1; //用来初始化感染者节点
int count1;
for(int i = 0; i < Network_Size; i++)
{
rate[i] = 0;
I_node[i] = 0;
kinf[i] = 0;
}
I = I0;
for(;;) //随机取初始感染者节点
{
if(I == 0)
{
break;
}
random_1 = genrand_int32()%Network_Size;
if(I_node[random_1] == 0)
{
I_node[random_1] = 1;
I--;
}
else
continue;
}
for(int i = 0; i < Network_Size; i++)
{
for(int j = 0; j < adjlist[i].size(); j++)
{
if(I_node[adjlist[i][j]] == 1)
{
kinf[i]++;
}
}
}
for(int i = 0; i < Network_Size; i++)
{
if(I_node[i] == 1)
{
rate[i] = g;
}
if(I_node[i] == 0)
{
rate[i] = r*kinf[i];
}
wt = wt + rate[i];
}
double v = genrand_real2();
dt = -log(1 - v)/wt;
I = I0;
count1 = 1;
X_vec_ind[0] = I;
t = dt;
for(;;)
{
double u = genrand_real2();
v = genrand_real2();
for(int i = 0; i < Network_Size; i++)
{
nk[i] = 0;
sum[i] = 0;
}
while((t > count1*interval)&&(count1 <= n_intervals))
{
X_vec_ind[count1] = I;
count1 ++;
}
if(t > tend)
{
break;
}
for(int i = 0; i < Network_Size; i++)
{
if(i == 0)
{
sum[i] = rate[i];
}
else
{
sum[i] = sum[i-1] + rate[i];
}
nk[i] = sum[i]/wt;
}
for(int i = 0; i < Network_Size; i++) //轮盘法选择节点发生状态转变
{
if(u < nk[0])
{
if(I_node[0] == 1)
{
I_node[0] = 0;
I--;
for(int j = 0; j < adjlist[i].size(); j++)
{
kinf[adjlist[0][j]]--;
}
break;
}
if(I_node[0] == 0)
{
I_node[0] = 1;
I++;
for(int j = 0; j < adjlist[0].size(); j++)
{
kinf[adjlist[0][j]]++;
}
break;
}
}
if(u > nk[i-1] && u < nk[i])
{
if(I_node[i] == 1)
{
I_node[i] = 0;
I--;
for(int j = 0; j < adjlist[i].size(); j++)
{
kinf[adjlist[i][j]]--;
}
break;
}
if(I_node[i] == 0)
{
I_node[i] = 1;
I++;
for(int j = 0; j < adjlist[i].size(); j++)
{
kinf[adjlist[i][j]]++;
}
break;
}
}
if(I == 0)
{
while (count1 <= n_intervals)
{
X_vec_ind[count1] = 0;
count1 ++;
}
break;
}
}
wt = 0;
for(int i = 0; i < Network_Size; i++)
{
if(I_node[i] == 1)
{
rate[i] = g;
}
if(I_node[i] == 0)
{
rate[i] = r*kinf[i];
}
wt = wt + rate[i];
}
dt = -log(1 - v)/wt;
t = t + dt;
}
X_vec.push_back(X_vec_ind);
}
double EXt;
for(int i=0; i <= n_intervals; i++) //取时间序列的平均值
{
EXt = 0;
for(int j=0; j < counter; j++)
EXt += X_vec[j][i];
EXt /= 1.0*counter*Network_Size;
outfile_ER_ext << EXt << "\n";
}
outfile_ER_ext<< endl;
}离散时间SIS模型与算法
而离散时间中，转变概率替代了转变速率，事件的转变时间由时间间隔 Δ𝑡 控制，其模型如下：
其中β'与μ'分别表示感染概率与恢复概率，值得注意的是，在大多数离散时间模型研究中，人们习惯使用β'=βΔ𝑡，μ'=μΔ𝑡，使用这种方法的一个弊端是连续时间与离散时间模型在参数较大时，两者无法映射(Physical Review E, 2016, 94(5):052125)。
离散时间模型算法通常采用同步更新方法：每个时间步骤中，所有个体都以同步方式更新其状态。状态为 S 的个体以概率 1-(1-β')^m 转变为状态 𝐼，而状态为 𝐼 的个体以概率 μ' 转变为状态 S 。整个过程迭代到系统达到稳定或系统中 𝐼 个体为零时停止。
C++算法实现
#include <iostream>
#include <ctime>
#include <stdlib.h>
#include <fstream>
#include <cmath>
#include "mt19937ar.c"
#include<algorithm>
#include <vector>
using namespace std;
int main()
{
int Network_Size = 10000;
double average_degree = 4;
double r = 0.2;
double g = 0.48;
double dt = 0.1;
double tend = 100;
double counter = 100;
int I0 = 100;
vector<int> degree(Network_Size);
unsigned long idum;
idum=(unsigned long)time(NULL);
init_genrand(idum);
vector<vector<int> > adjlist(Network_Size);
ofstream outfile_ER_S_ext;
double r1;
outfile_ER_S_ext.open("ER_S_ext.txt");
double Sum_degree = 0;
int edge = 0.5*Network_Size*average_degree;
for(;;)
{
double judge = 0;
int random_1 = genrand_int32()%Network_Size;
int random_2 = genrand_int32()%Network_Size;
if(random_1 != random_2)
{
for(int i = 0; i < adjlist[random_1].size(); i++)
{
if(adjlist[random_1][i] == random_2)
{
judge = 1;
break;
}
}
if(judge == 1)
{
continue;
}
else
{
adjlist[random_1].push_back(random_2);
degree[random_1]++;
adjlist[random_2].push_back(random_1);
degree[random_2]++;
edge--;
}
}
if(edge == 0)
break;
}
double rs;
double gs;
double t = 0;
rs = r*dt;
gs = g*dt;
int I;
int random_1; //用来初始化感染者节点
vector<int> kinf(Network_Size,0);
vector<int> I_node(Network_Size,0);
vector<int> kinf_old(Network_Size);
int n_intervals = tend/dt;
double interval = tend/(1.0*n_intervals);
vector<int> X_vec_ind(n_intervals + 1);
vector<vector<int> > X_vec;
int count1;
for(int n = 0; n < counter; n++)
{
for(int i = 0; i < Network_Size; i++)
{
kinf[i] = 0;
I_node[i] = 0;
}
I = I0;
for(;;) //随机取初始感染者节点
{
if(I == 0)
{
break;
}
random_1 = genrand_int32()%Network_Size;
if(degree[random_1] == 0)
{
continue;
}
if(I_node[random_1] == 0)
{
I_node[random_1] = 1;
I--;
}
else
continue;
}
for(int i = 0; i < Network_Size; i++)
{
for(int j = 0; j < adjlist[i].size(); j++)
{
if(I_node[adjlist[i][j]] == 1)
{
kinf[i]++;
}
}
}
t = 0;
I = I0;
X_vec_ind[0] = I;
count1 = 1;
for(;;)
{
for(int i = 0; i < Network_Size; i++)
{
kinf_old[i] = kinf[i];
}
t = t + dt;
while((t > count1*interval)&&(count1 <= n_intervals))
{
X_vec_ind[count1] = I;
count1 ++;
}
if(t > tend)
break;
if(I == 0)
{
while (count1 <= n_intervals)
{
X_vec_ind[count1] = 0;
count1 ++;
}
break;
}
double u;
for(int i = 0; i < Network_Size; i++)
{
u = genrand_real2();
if(I_node[i] == 1)
{
if(u < gs)
{
I_node[i] = 0;
I--;
for(int j = 0; j < adjlist[i].size(); j++)
{
kinf[adjlist[i][j]]--;
}
}
}
else if(I_node[i] == 0)
{
if(u < 1 - pow((1 - rs),kinf_old[i]))
{
I_node[i] = 1;
I++;
for(int j = 0; j < adjlist[i].size(); j++)
{
kinf[adjlist[i][j]]++;
}
}
}
}
}
X_vec.push_back(X_vec_ind);
}
double EXt;
for(int i=0; i <= n_intervals; i++)
{
EXt = 0;
for(int j=0; j < counter; j++)
EXt = EXt + X_vec[j][i];
EXt = EXt/(counter*Network_Size);
outfile_ER_S_ext << EXt << "\n";
}
outfile_ER_S_ext << endl;
return 0;
}
结语
对于各种网络的实现方法大同小异，但对于不同的网络需要单独优化程序，例如在全连接网络中，每个节点的感染者邻居数是相同的，因此可以极大的简化抽样的步骤。上面这两个程序是我刚入门时参考论文中大佬给出的程序写的，在后面的独立研究中也很少使用到，如有错误，请指正，谢谢。
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
