python2.7 SIRSシミュレーションをアニメーションを使わずに実装したい

python初心者です。

###前提・実現したいこと
ネットワーク攻撃シミュレーション実験のために、SIRS感染モデルをpython上で実装し、シミュレーションしようとしています。
verはpyrhon2.7です。
当初はgraph_toolを利用してアニメーションを用いて実装していましたが、データを取る時に一々ウィンドウを閉じなければならなかったので、アニメーションを使わずに実装することにしました。
以下に、ソースコードを示します。

###該当のソースコード

python
1#! /usr/local/bin/env python
2# coding:utf-8
3
4from numpy.random import *
5from graph_tool.all import *
6import networkx as nx
7import sys, os, os.path
8
9import create_graph as g # create_graphはネットワークグラフを作るためのモジュール
10G=g.mgraph() # グラフを保持する場所
11# グラフのノード数は１００
12
13#グラフの読み込みが終わったので、ここからはグラフのノードについて操作
14g2=nx.Graph()
15for e in G.edges():
16    g2.add_edge(int(e.source()), int(e.target()))
17
18removed = G.new_vertex_property("bool")
19
20beta = 0.1  # 感染率
21s = 0.01    # I->Sになる確率(回復率)
22r = 0.01    # I->Rになる確率(死亡率)
23
24S = [1, 1, 1, 1]          # 状態S  
25I = [0, 0, 0, 0]          # 状態I
26R = [0.8, 0.8, 0.8, 0.8]  # 状態R
27
28
29# 最初は状態Sで初期化
30state = G.new_vertex_property("vector<double>")
31for v in G.vertices():
32    state[v] = S
33
34# ランダムに10個のノードを感染
35vs = list(G.vertices())
36for var in range(0,10):
37    i = vs[randint(0, len(vs))]
38    if(state[i] != I):
39        state[i] = I
40        
41newly_infected = G.new_vertex_property("bool") # 新たに感染したノード
42
43# 状態遷移させる関数
44def update_state():
45    newly_infected.a = False
46    removed.a = False
47
48    # visit the nodes in random order
49    vs = list(G.vertices())
50    shuffle(vs)
51    scount=0      # 状態Sを数える
52    icount=0      # 状態Iを数える
53    rcount=0      # 状態Rを数える
54    for v in vs:
55        if state[v] == I:
56            icount += 1
57            if random() < r:
58                state[v] = R
59                try:
60                    g2.remove_node(v)
61                except:
62                    pass
63            elif random() < s :
64                state[v] = S
65        elif state[v] == S:
66            scount+=1
67            for w in v.out_neighbours():
68                if state[w] == I:
69                    if (random()<beta):
70                        state[v] = I
71                        newly_infected[v] = True
72                        break
73        if state[v] == R:
74            rcount += 1
75            removed[v] = True
76        
77        if icount==0:
78            GC = sorted(nx.connected_component_subgraphs(g2), key=len, reverse=True)
79            GC=len(GC[0])
80            print ("s:",scount)
81            print ("i:",icount)
82            print ("r:",rcount)
83
84
85update_state()

このソースコードを動かすと、特定のエラーメッセージは出なかったのですが、思った通りの挙動になりませんでした。
想定では以下のように、

('s:', 90)
('i:', 10)
('r:', 0)
('s:', 84)
('i:', 14)
('r:', 2)
・・・
・・・
('s:', 76)
('i:', 0)
('r:', 24)

というような遷移をします(実際に、アニメーションを実装した状態では上記のような結果が出力されました。)。
ところが、先ほどのソースコードを動かすと以下のような出力になりました。

('s:', 1)
('i:', 0)
('r:', 0)
('s:', 2)
('i:', 0)
('r:', 0)
・・・
・・・
('s:', 32)
('i:', 0)
('r:', 0)

この出力のおかしな点は3点あります。
・状態遷移させる前に、ランダムに１０個感染させているのに I の数が0個になっている。
・確率で遷移させているのに I と R の数が増えず、S の数が1個ずつ増えている(ように見える)。
・そもそも上記のような結果が出ない時がある(実行される時とされない時がある)。

###試したこと
関数の引数の問題かと思い、引数を指定してみましたが結果は変わりませんでした。

アニメーション部分のコードを丸ごとコメントアウトした以外は特に書き換えた部分はなく、なぜ思ったような挙動にならないかわからず困っています。
どうすれば思い通りの挙動になるでしょうか？