matplotlibで複数のグラフの平均をとってグラフを描画する方法がわかりません。

Question

実現したいこと

• 複数のグラフの平均をとって一つのグラフを出力する。

前提

ここに質問の内容を詳しく書いてください。

pythonで以下のようにネットワーク上のSIRモデルをシミュレーションすることを目指しております。最後、10回シミュレーションを行ってデータを重ねてグラフを描画するところまでは完了しました。しかし、その10個のグラフの平均をとって一つのグラフを作成する方法がよくわかりません。

発生している問題・エラーメッセージ

白紙のグラフが出力されました。

該当のソースコード

import matplotlib.pyplot as plt
import numpy as np
import networkx as nx

#ネットワークの設定。今回はBAモデル
N=10000 #人口
M=5 #次数
G=nx.barabasi_albert_graph(N,M)

#SIRモデルのシミュレーションを定義
class SIR:
    def __init__(self, G, beta, gamma): #ネットワーク、感染率、回復率を定義
        self.G = G
        self.beta = beta
        self.gamma = gamma
        
    def run(self, seed=[], num_steps = 1, sentinels = [], immunization_rate = 0.0, immunized_nodes = []):
        #ネットワーク上の初期状態での免疫保持者をランダムに決定する。今回は0人。 
        if len(immunized_nodes) == 0:
            immunized = set(np.random.choice(self.G.nodes(), 
                                             size=int(immunization_rate*len(self.G.nodes())), 
                                             replace=False))
        else:
            immunized = immunized_nodes

        # 初期状態でのIをランダムで決定する。
        if len(seed) == 0:
            seed = [np.random.choice(list(set(self.G.nodes()).difference(immunized)))]
        
        #S,I,Rの初期状態
        I_set = set(seed)
        S_set = set(self.G.nodes()).difference(I_set).difference(immunized)
        R_set = set()
        
        t = 0
        
        StoI = set(seed)
        ItoR = set()
        
        sentinels_t = {}
        for sen in sentinels:
            sentinels_t[sen] = 0
        
        #I>0の場合、感染症が伝播
        while len(I_set) > 0:
            for i in I_set.copy():
                for s in set(self.G.neighbors(i)).intersection(S_set).copy():
                    if np.random.uniform() < self.beta:
                        S_set.remove(s)
                        I_set.add(s)
                        StoI.add(s)
                        
                        if sentinels_t.get(s) != None:
                            sentinels_t[s] = t
                        
                #回復率に従ってIをRに変化
                if np.random.uniform() < self.gamma:
                    I_set.remove(i)
                    R_set.add(i)
                    ItoR.add(i)
    
            t += 1
            if t % num_steps == 0 or len(I_set) == 0:
                yield({'t': t, 'S':S_set, 'I':I_set, 'R':R_set, 'StoI':StoI, 'ItoR':ItoR, 'sentinels': sentinels_t})
beta = 0.2 #感染率
gamma = 0.1 #回復率

peak_times = []

#複数のデータを一つのグラフにまとめる。
plt.figure()
ax = plt.subplot(111)
ax.set_prop_cycle('color', ['green', 'red', 'blue'])

#試行回数を決定。L回シミュレーションしてデータを重ねる。
L=10
for i in range(0,L):
    sir = SIR(G, beta = beta, gamma = gamma)
    ave = []
    sum=0
    for r in sir.run(num_steps=1):
        sum_S=sum+len(r['S'])
        sum_I=sum+len(r['I'])
        sum_R=sum+len(r['R'])
        ave.append([sum_S, sum_I, sum_R])
        

    #データをプロット
    plt.plot(ave)
    
#結果をプロット
plt.title("Infection On Network")
plt.legend(['S', 'I', 'R'])
plt.xlabel("day")
plt.ylabel("Population")
plt.show()

試したこと

numpyで平均値を取る関数を用いたりしましたが、うまくいきませんでした。
有識者の方、どうか解決法をご教授いただけませんでしょうか。

Comments

[8524ba23]

まずは１回だけシミュレーションした結果は正しく描画できるでしょうか？
提示コードをみる限り
- sum=0 はつねに0のままのようだがそれでよい？
-「sum_S=～」行のインデントは正しいか？
- resは利用されていないがそれでよい？
など、ところどころ意図が読めない部分がありますので、そのあたりの説明を追記すると回答得られやすくなるかと思います。

[onomatobee]

はい、1回だけ行った場合は正しく描画できました。
10回ほどシミュレーションを行った後、その平均を取る方法をご教授いただければと思います。

[onomatobee]

ご質問ありがとうございます。
私なりに当該箇所を修正させていただきました。

[melian]

sum_S=sum+len(r['S'])
sum_I=sum+len(r['I'])
sum_R=sum+len(r['R'])

の部分ですが、これは、以下の様にするのではないでしょうか？

sum_S = sum_I = sum_R = 0
for r in sir.run(num_steps=1):
　sum_S += len(r['S'])
　sum_I += len(r['I'])
　sum_R += len(r['R'])
total = sum_S + sum_I + sum_R
ave.append([sum_S/total, sum_I/total, sum_R/total])

[onomatobee]

ありがとうございます。
ぜひ試させていただきたいと思います。

[onomatobee]

試してみたのですが、どうやら時間tに関しても平均値が取られてしまっている様です‥。

Answer 1

#試行回数を決定。
L = 10
sumup, length = [], []
for _ in range(L):
    sir = SIR(G, beta = beta, gamma = gamma)
    for i, r in enumerate(sir.run(num_steps=1)):
        sum_S, sum_I, sum_R = len(r['S']), len(r['I']), len(r['R'])
        if i < len(sumup):
            sumup[i][0] += sum_S
            sumup[i][1] += sum_I
            sumup[i][2] += sum_R
            length[i] += 1
        else:
            sumup.append([sum_S, sum_I, sum_R])
            length.append(1)

#データをプロット
ave = [[v / l for v in s] for s, l in zip(sumup, length)]
plt.plot(ave)

Answer 2

ご回答本当にありがとうございます。
早速試させていただいたのですが、やはり10個のデータが別々に出力されてしまう様です…。

Answer 3

10個のグラフの平均をとって一つのグラフを作成する方法

numpyやmatplotlibに便利な関数は存在しないとおもうので、地道に各試行結果を計算するほかないと思います。

７５行目以下を下記コードに変更するのでいかがでしょうか？
※平均のグラフを出すことを目標としたので、グラフの色などの設定は適当です。

import pandas as pd
L=10
res_for_ave = [] #各試行の結果を保存しておくためのリスト
for i in range(0,L):
    sir = SIR(G, beta = beta, gamma = gamma)
    res = []
    ave = []
    sum=0
    for r in sir.run(num_steps=1):
        res.append([len(r['S']), len(r['I']), len(r['R'])])
    plt.plot(res)
    res_for_ave.append(res) #試行結果を保存

# pandasで保存したデータを加工し、平均を算出
df = pd.DataFrame(res_for_ave)
df = df.T
df = df.fillna(0)
x = np.arange(0,len(df)) #時間
ave_s = [] #S
ave_i =[] #I
ave_r = [] #R
for row in range(len(df)): #時間でLoop
    count_sum_s = 0
    count_sum_i = 0
    count_sum_r = 0
    count_len = L
    for col in df.columns: #試行回数分でLoop
        if df[col][row] == 0:
            count_len -= 1
        else:
            count_sum_s += df[col][row][0] #S
            count_sum_i += df[col][row][1] #I
            count_sum_r += df[col][row][2] #r
    ave_s.append(count_sum_s/count_len)
    ave_i.append(count_sum_i/count_len)
    ave_r.append(count_sum_r/count_len)
plt.plot(x,ave_s, label = "S")
plt.plot(x,ave_i, label = "I")
plt.plot(x,ave_r, label = "R")
plt.title("Infection On Network")
plt.legend(['S', 'I', 'R'])
plt.xlabel("day")
plt.ylabel("Population")
plt.show()

Comments

[onomatobee]

再度試したところ、無事解決することができました！
本当にありがとうございます！