吸収マルコフ連鎖の紹介

1
2
from scipy.stats import chisquare
chisquare([69, 53, 14, 6, 26], f_exp=v_7)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
import numpy as np

class Revere_Plan_Sumulator:

    def __init__(self, population_size, initial_news_holder_ratio, maximum_state, eps = 1e-18):
        """ 
        Args
            population_size
                コミュニティの構成人数
                そのコミュニティを対象にニュースが広まる程度を観察する

            initial_news_holder_ratio
                ニュースを初めに伝達した人数の割合＝伝達人数割合初期値

            maximum_state
                コニュニティの構成員のstateの種類の上限値
                (例) もし、maximum_state = 3 ならば、コニュニティの構成員のstateは三人以上を挟んでニュースを又聞きしたまでとなる

        Outputs
            prob_vector
                状態割合変数
                S_0, S_1, ...., S_{maximum_state}, G_0, G_1, ...., G_{maximum_state}, I

        How to use it
            EXAMPLE: 
                newsの各期の新規伝達社人数の推移のplot

            CODE
                community_size = 210
                initial_news_holder = 0.2
                Maximum_state = 4
                iter_number = 1000

                TestClass = Revere_Plan_Sumulator(community_size, initial_news_holder, Maximum_state)
                prob_vector, spread_sequence = TestClass.do_simulation(iter_number)
                
                plt.plot(spread_sequence[1:]) 
        """

        ## input
        self.population_size = population_size
        self.initial_news_holder = initial_news_holder_ratio
        self.maximum_state= maximum_state
        self.eps = eps
        self.news_holder_sequence = None

        # create an initial vector
        field_number = 2 * (self.maximum_state + 1) + 1
        self.prob_vector = np.zeros(field_number)
        self.prob_vector[maximum_state+1], self.prob_vector[-1] = initial_news_holder_ratio, 1 - initial_news_holder_ratio
        
        # create an transition matrix
        self.transition_matrix = np.zeros((field_number, field_number))
        self.transition_matrix[:self.maximum_state + 1, :self.maximum_state + 1] = np.eye(self.maximum_state + 1)

    def update_transition_matrix(self):
        """
        Description
            今季の遷移行列の計算
        
        Note
            1. (吸収状態, 吸収状態)の遷移は単位行列
            2.  (吸収状態, 非吸収状態)の遷移は0
            3. Q: (非吸収状態, 非吸収状態)の遷移は基本的には１期前の「Naive」人数割合が対角要素となる行列
            4.  R: (非吸収状態, 吸収状態)の遷移は「単位行列 - (非吸収状態, 非吸収状態)の行列」
            5. last_row: 「naive」のspreaderへの遷移は１期前の「spreader」人数割合をずらして作成
        
        REMARKS
            (非吸収状態, 非吸収状態)の遷移はノートと照らし合わせると、１期前にG_iが0以上ならば 
            G_i→G_i の遷移が１期前の「Naive」人数割合となり、それ以外ならば 0 とすべきだが
            計算上問題ない(１期前のベクトル要素が0ならば何があっても0)なので、実装からは
            除去
         
            which_element_update_vector = np.where(self.prob_vector[self.maximum_state+1:-1] > self.eps, 1, 0) ##ここは計算上なくても問題ない
            Q = np.eye(self.maximum_state+1) * (which_element_update_vector * self.prob_vector[-1]) ## convert to spreader
            R = np.eye(self.maximum_state+1) * (which_element_update_vector * (1 - self.prob_vector[-1])) ## convert to suppressor
        """

        Q = np.eye(self.maximum_state+1) * self.prob_vector[-1] ## convert to spreader
        R = np.eye(self.maximum_state+1) * (1 - self.prob_vector[-1]) ## convert to suppressor
        
        last_row = self.prob_vector[self.maximum_state+1:-1].copy() ## 配列コピー
        last_row[-2] = last_row[-2] + last_row[-1]
        last_row[-1] = max(1 - sum(last_row[:-1]),0)

        self.transition_matrix[self.maximum_state+1:-1, :self.maximum_state+1] = R
        self.transition_matrix[self.maximum_state+1:-1, self.maximum_state+1:-1] = Q
        self.transition_matrix[-1,self.maximum_state+2:] = last_row

    def update_prob_vector(self):
        """  
        Description
            今季の状態割合ベクトル
            = １期前の状態割合ベクトル × 確率遷移行列

        """
        ## 今季の確率遷移行列のupdate
        self.update_transition_matrix()
        
        
        return np.matmul(self.prob_vector, self.transition_matrix)

    def do_simulation(self, iter_limit):
        """  
        Args
            iter_limit
                何回目の出会いまで計算するかをsimulation
        
        """
        count_iter = 0
        
        self.news_holder_sequence = np.zeros(iter_limit+1)
        self.news_holder_sequence[0] = self.initial_news_holder

        while count_iter < iter_limit:
            count_iter += 1
            prob_vec = self.update_prob_vector()
            
            news_holder_increment = self.prob_vector[-1] - prob_vec[-1]
            self.news_holder_sequence[count_iter] = news_holder_increment
            
            if abs(news_holder_increment) < self.eps:
                break
            else:
                self.prob_vector = prob_vec

        return self.prob_vector, self.news_holder_sequence[:count_iter]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chisquare

plt.rcParams['font.family'] = 'Meiryo'

M = 4
TestClass = Revere_Plan_Sumulator(210, 0.2, M)
prob_vector, spread_sequence = TestClass.do_simulation(1000)

plt.figure(figsize=(10, 8))
plt.plot(spread_sequence[1:])
plt.title('一過性のNews流行の変化')
plt.xlabel('Markov chainの各期間')
plt.ylabel('新規News伝聞者割合');