用隐马尔可夫模型预测股票（大雾）

 1    def __init__(self, trans_pbb: np.ndarray, emit_pbb: np.ndarray, init_pbb: np.ndarray):
 2        """初始化隐马尔可夫模型
 3
 4        Args:
 5            trans_pbb (np.ndarray): 转移概率矩阵
 6            emit_pbb (np.ndarray): 输出概率矩阵
 7            init_pbb (np.ndarray): 初始概率分布
 8        """
 9        self.trans_pbb = trans_pbb
10        self.emit_pbb = emit_pbb
11        self.init_pbb = init_pbb
12        self.n_states = len(init_pbb)

 1    n_seq = 28 # 预测长度
 2    n_turn = 100 # 随机次数
 3    
 4    # == 使用前向算法，估计最后一个状态的概率分布 ==
 5    # 由于已经在前面调用过，所以我们直接使用 alpha 计算结果
 6    alpha_last = alpha[:, 99]
 7    output_real = ob_seq[100:128]  # 实际值
 8
 9    # 最后状态。根据前向算法得出
10    forward_last_state = np.argmax(alpha_last)
11    price_seq_real = _reduce_output(output_real) # 真实价格序列
12    acc_seq = np.zeros(n_turn)
13    # 进行 n_turn 次模拟
14    for turn in range(n_turn):
15        last_state = forward_last_state
16        output_fit = np.zeros(n_seq)  # 预测值
17        for t in range(n_seq):
18            # 按照状态转移概率分布，抽取一个状态
19            state = np.random.choice(
20                np.arange(3), p=model.trans_pbb[last_state])
21            # 按照此状态到各个观测值的概率，抽取一个观测值
22            output = np.random.choice(np.arange(5), p=model.emit_pbb[state, :])
23            output_fit[t] = output
24
25            last_state = state
26        from sklearn.metrics import precision_score
27        acc_seq[turn] = precision_score(output_real, output_fit, average='micro')
28        
29
30    # 输出拟合度的均值和标准方差
31    print("acc_seq=")
32    print(acc_seq)
33    print("acc_seq.mean=")
34    print(acc_seq.mean())
35    print("acc_seq.std=")
36    print(acc_seq.std())

 1
 2```python
 3    def _reduce_output(output_list: np.ndarray):
 4        """根据输出的状态序列，从初始值开始累加，形成一个（拟）价格序列
 5
 6        Args:            
 7            output_list (np.ndarray): 输出序列
 8
 9        Returns:
10            ndarray: （拟）价格序列
11        """
12        init_val = 1
13        # 从状态映射到涨跌幅。百分比。
14        price_map = {
15            0: -2,  # large drop
16            1: -1.5,  # drop
17            2: 0,  # no change
18            3: 1.5,  # small rise
19            4: 2  # large rise
20        }
21        enable_price_map = True
22        ret = np.zeros(len(output_list))
23        ret[0] = init_val
24        if enable_price_map:
25            for i in range(1, len(output_list)):
26                ret[i] = ret[i - 1] * (1 + price_map[output_list[i]] / 100)
27        else:
28            for i in range(1, len(output_list)):
29                ret[i] = ret[i - 1] + output_list[i] - 2
30        return ret

1from scipy.io import loadmat
2mat = loadmat(file_mat, squeeze_me=True)

 1    # trans_bb[qi][qj] 表示状态 qi 转移到 qj 的概率
 2    trans_bb = np.array(
 3        [[0.26596786, 0.64534235, 0.08868979],
 4         [0.49997787, 0.49788423, 0.00213789],
 5         [0.02295626, 0.00247515, 0.97456859]]
 6    )
 7
 8    # emit_bb[qi][oj] 表示状态 qi 下观测到 oj 的概率
 9    emit_bb = np.array(
10        [[4.37786679e-03, 3.56626266e-02, 9.04597402e-01, 5.18816277e-02, 3.48047712e-03],
11         [1.17931200e-01, 2.36821594e-01, 2.51012581e-01,
12             2.21650346e-01, 1.72584279e-01],
13         [9.45849521e-03, 8.32174208e-02, 8.54638596e-01, 5.26834471e-02, 2.04043929e-06]]
14    )
15
16    # prior_bb[qi] 表示状态 qi 的先验概率。即初始取 qi 的概率。
17    init_pbb = np.array([2.52355678e-27, 7.86802271e-01, 2.13197729e-01])
18
19    # ob_seq[t] 表示第 t 个时刻的观测值
20    ob_seq = np.array([2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3,
21                       3, 3, 3, 3, 3, 4, 3, 5, 4, 3, 3, 3, 2, 3, 4, 3, 2, 3, 4, 3, 3, 4,
22                       3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 2, 3, 3, 3, 3, 3, 3, 2, 3, 3, 3,
23                       3, 3, 3, 3, 2, 3, 3, 3, 4, 2, 5, 3, 3, 2, 3, 3, 3, 3, 4, 3, 3, 3,
24                       3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 3,
25                       4, 3, 3, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 4, 3, 3, 3], dtype=np.uint8)
26
27    # 因为索引从 0 开始，所以修正一下偏移量
28    ob_seq = ob_seq - 1

 1    def calc_forward_pbb(self, obs_seq: int):
 2        """计算前向概率矩阵
 3
 4        Args:
 5            obs_seq (int): 观测序列
 6
 7        Returns:
 8            np.ndarrya: 前向概率矩阵，维度为 (n_states, len(obs_seq))
 9        """
10        n_obs = len(obs_seq)
11        forward_pbb_mat = np.zeros((self.n_states, n_obs))
12        forward_pbb_mat[:, 0] = self.init_pbb * self.emit_pbb[:, 0]
13        for t in range(1, n_obs):
14            for i in range(self.n_states):
15                forward_pbb_mat[i, t] = np.dot(
16                    forward_pbb_mat[:, t - 1], self.trans_pbb[:, i]
17                ) * self.emit_pbb[i, obs_seq[t]]
18
19        return forward_pbb_mat

1forward_pbb_mat[i, t] = 
2  np.dot(forward_pbb_mat[:, t - 1], self.trans_pbb[:, i]) # 表示时间前进一步之后产生状态 i 的概率
3                          * self.emit_pbb[i, obs_seq[t]]  # 表示状态 i 产生 t 对应的观测值的概率
4                                                          # 二者相乘，表示时间前进一步之后，状态 i 产生观测值 obs_seq[t] 的概率

 1    def calc_backward_pbb(self, ob_seq: list[int]):
 2        """计算后向概率矩阵
 3
 4        Args:
 5            ob_seq (list[int]): 观测序列
 6        """
 7        n_obs = len(ob_seq)
 8        backward_pbb_mat = np.zeros((self.n_states, n_obs))
 9        backward_pbb_mat[:, -1] = 1
10
11        for t in reversed(range(n_obs - 1)):
12            for i in range(self.n_states):
13                backward_pbb_mat[i, t] = np.sum(
14                    backward_pbb_mat[:, t + 1] * self.trans_pbb[i, :] *
15                    self.emit_pbb[:, ob_seq[t+1]]
16                )
17        return backward_pbb_mat

1backward_pbb_mat[i, t] = np.sum(
2    backward_pbb_mat[:, t + 1] * # 表示下一时刻的各状态后向概率
3    self.trans_pbb[i, :] *       # 表示从 i 状态转移到其它状态的概率
4    self.emit_pbb[:, ob_seq[t+1]]# 表示发射出 t+1 时输出的概率
5)

 1    def viterbi_decode(self, ob_seq: list[int]):
 2        """计算隐藏状态序列（维特比路径）
 3
 4        Args:
 5            ob_seq (list[int]): 观测序列
 6        """
 7        n_obs = len(ob_seq)
 8        prev = np.zeros((n_obs - 1, self.n_states))
 9        
10        # 构造 DP 数组
11        # dp[i, t] 表示在 t 时刻，到达 i 状态的观测序列的最大概率
12        dp = np.zeros((self.n_states, n_obs))
13        dp[:, 0] = self.init_pbb * self.emit_pbb[:, ob_seq[0]]
14        for t in range(1, n_obs):
15            for i in range(self.n_states):
16                seq_pbb = dp[:, t-1] * self.trans_pbb[:, i] * \
17                    self.emit_pbb[i, ob_seq[t]]
18                prev[t-1, i] = np.argmax(seq_pbb)
19                dp[i, t] = np.max(seq_pbb)
20
21        # 构造维特比路径
22        def build_viterbi_path(prev, last_state):
23            T = len(prev)
24            yield(last_state)
25            for i in range(T-1, -1, -1):
26                v = int(prev[i, last_state])
27                yield(v)
28                last_state = v
29                
30        viterbi_path = np.array(list(build_viterbi_path(prev, np.argmax(dp[:, -1]))))
31        return viterbi_path

  1from matplotlib import pyplot as plt
  2import numpy as np
  3
  4
  5class HMM():
  6
  7    def __init__(self, trans_pbb: np.ndarray, emit_pbb: np.ndarray, init_pbb: np.ndarray):
  8        """初始化隐马尔可夫模型
  9
 10        Args:
 11            trans_pbb (np.ndarray): 转移概率矩阵
 12            emit_pbb (np.ndarray): 输出概率矩阵
 13            init_pbb (np.ndarray): 初始概率分布
 14        """
 15        self.trans_pbb = trans_pbb
 16        self.emit_pbb = emit_pbb
 17        self.init_pbb = init_pbb
 18        self.n_states = len(init_pbb)
 19
 20    def calc_forward_pbb(self, ob_seq:  list[int]):
 21        """计算前向概率矩阵
 22
 23        forward_pbb_mat[i, t] 表示隐藏态为 i 条件下，在 0~t 时刻的前向概率。也即，经过 0~t 步，生成观测序列 ob_seq ，并最终落在 i 状态的概率。
 24
 25        Args:
 26            ob_seq (list[int]): 观测序列
 27
 28        Returns:
 29            np.ndarrya: 前向概率矩阵，维度为 (n_states, len(ob_seq))
 30        """
 31        n_obs = len(ob_seq)
 32        forward_pbb_mat = np.zeros((self.n_states, n_obs))
 33        forward_pbb_mat[:, 0] = self.init_pbb * self.emit_pbb[:, 0]
 34        for t in range(1, n_obs):
 35            for i in range(self.n_states):
 36                forward_pbb_mat[i, t] = np.dot(
 37                    forward_pbb_mat[:, t - 1], self.trans_pbb[:, i]
 38                ) * self.emit_pbb[i, ob_seq[t]]
 39
 40        return forward_pbb_mat
 41
 42    def calc_backward_pbb(self, ob_seq: list[int]):
 43        """计算后向概率矩阵
 44
 45        backward_pbb_mat[i, t] 表示隐藏态为 i 条件下，在 t~T 时刻的后向概率。也即，从 i 状态开始，经过 t~T 步，生成观测序列 ob_seq 的概率。
 46
 47        Args:
 48            ob_seq (list[int]): 观测序列
 49        """
 50        n_obs = len(ob_seq)
 51        backward_pbb_mat = np.zeros((self.n_states, n_obs))
 52        backward_pbb_mat[:, -1] = 1
 53
 54        for t in reversed(range(n_obs - 1)):
 55            for i in range(self.n_states):
 56                backward_pbb_mat[i, t] = np.sum(
 57                    backward_pbb_mat[:, t + 1] * self.trans_pbb[i, :] *
 58                    self.emit_pbb[:, ob_seq[t+1]]
 59                )
 60        return backward_pbb_mat
 61
 62    def viterbi_decode(self, ob_seq: list[int]):
 63        """计算隐藏状态序列（维特比路径）
 64
 65        Args:
 66            ob_seq (list[int]): 观测序列
 67        """
 68        n_obs = len(ob_seq)
 69        prev = np.zeros((n_obs - 1, self.n_states))
 70
 71        # 构造 DP 数组
 72        # dp[i, t] 表示在 t 时刻，到达 i 状态的观测序列的最大概率
 73        dp = np.zeros((self.n_states, n_obs))
 74        dp[:, 0] = self.init_pbb * self.emit_pbb[:, ob_seq[0]]
 75        for t in range(1, n_obs):
 76            for i in range(self.n_states):
 77                seq_pbb = dp[:, t-1] * self.trans_pbb[:, i] * \
 78                    self.emit_pbb[i, ob_seq[t]]
 79                prev[t-1, i] = np.argmax(seq_pbb)
 80                dp[i, t] = np.max(seq_pbb)
 81
 82        # 构造维特比路径
 83        def build_viterbi_path(prev, last_state):
 84            T = len(prev)
 85            yield(last_state)
 86            for i in range(T-1, -1, -1):
 87                v = int(prev[i, last_state])
 88                yield(v)
 89                last_state = v
 90
 91        viterbi_path = np.array(
 92            list(build_viterbi_path(prev, np.argmax(dp[:, -1]))))
 93        return viterbi_path
 94
 95
 96def draw_state_trans(state_seq, fname):
 97    seq_time = np.arange(len(state_seq))
 98    fg, ax = plt.subplots(figsize=(10, 2))
 99    ax.set(xlabel='time', ylabel='state', title='Hidden states')
100    state0 = np.ma.masked_where(state_seq == 0, state_seq)
101    state1 = np.ma.masked_where(state_seq == 1, state_seq)
102    state2 = np.ma.masked_where(state_seq == 2, state_seq)
103    ax.step(seq_time, state0, '-', label='state 0 (bear)')
104    ax.step(seq_time, state1, '-', label='state 1 (stable)')
105    ax.step(seq_time, state2, '-', label='state 2 (bull)')
106    ax.legend()
107    plt.savefig(fname)
108
109def draw_pbb_diagram(gamma, fname):
110    seq_time = np.arange(len(gamma[0]))
111    fg, ax = plt.subplots(figsize=(10, 2))
112    ax.set(xlabel='time', ylabel='probability', title='Probability Dist.')
113    state0data = gamma[0]
114    state1data = gamma[1]
115    state2data = gamma[2]
116    ax.plot(seq_time, state0data, '-', label='state 0 (bear)')
117    ax.plot(seq_time, state1data, '-', label='state 1 (stable)')
118    ax.plot(seq_time, state2data, '-', label='state 2 (bull)')
119    ax.legend()
120    plt.savefig(fname)
121
122
123def main():
124
125    # trans_bb[qi][qj] 表示状态 qi 转移到 qj 的概率
126    trans_bb = np.array(
127        [[0.26596786, 0.64534235, 0.08868979],
128         [0.49997787, 0.49788423, 0.00213789],
129         [0.02295626, 0.00247515, 0.97456859]]
130    )
131
132    # emit_bb[qi][oj] 表示状态 qi 下观测到 oj 的概率
133    emit_bb = np.array(
134        [[4.37786679e-03, 3.56626266e-02, 9.04597402e-01, 5.18816277e-02, 3.48047712e-03],
135         [1.17931200e-01, 2.36821594e-01, 2.51012581e-01,
136             2.21650346e-01, 1.72584279e-01],
137         [9.45849521e-03, 8.32174208e-02, 8.54638596e-01, 5.26834471e-02, 2.04043929e-06]]
138    )
139
140    # prior_bb[qi] 表示状态 qi 的先验概率。即初始取 qi 的概率。
141    init_pbb = np.array([2.52355678e-27, 7.86802271e-01, 2.13197729e-01])
142
143    # ob_seq[t] 表示第 t 个时刻的观测值
144    ob_seq = np.array([2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3,
145                       3, 3, 3, 3, 3, 4, 3, 5, 4, 3, 3, 3, 2, 3, 4, 3, 2, 3, 4, 3, 3, 4,
146                       3, 3, 3, 3, 3, 3, 3, 3, 2, 3, 3, 2, 3, 3, 3, 3, 3, 3, 2, 3, 3, 3,
147                       3, 3, 3, 3, 2, 3, 3, 3, 4, 2, 5, 3, 3, 2, 3, 3, 3, 3, 4, 3, 3, 3,
148                       3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 3,
149                       4, 3, 3, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 4, 3, 3, 3], dtype=np.uint8)
150
151    # 因为索引从 0 开始，所以修正一下偏移量
152    ob_seq = ob_seq - 1
153
154    model = HMM(trans_bb, emit_bb, init_pbb)
155    ob_seq_100 = ob_seq[:100]
156    alpha = (model.calc_forward_pbb(ob_seq_100))
157    beta = (model.calc_backward_pbb(ob_seq_100))
158    print("alpha=")
159    print(alpha)
160    print("beta=")
161    print(beta)
162
163    # 通过求各时刻个状态的最大值，求状态序列
164    seq_gamma = alpha * beta
165    seq_gamma_argmax = np.argmax(seq_gamma, axis=0)
166    # 绘制
167    print("gamma=")
168    print(seq_gamma)
169    print("gamma(argmax)=")
170    print(seq_gamma_argmax)
171
172    draw_state_trans(seq_gamma_argmax, fname='problem_1.png')
173    draw_pbb_diagram(alpha * beta, fname='problem_1_2.png')
174
175    viterbi_seq = model.viterbi_decode(ob_seq_100)
176    print("viterbi_seq=")
177    print(viterbi_seq)
178    draw_state_trans(viterbi_seq, fname='problem_2.png')
179
180    def _reduce_output(output_list: np.ndarray):
181        """根据输出的状态序列，从初始值开始累加，形成一个（拟）价格序列
182
183        Args:            
184            output_list (np.ndarray): 输出序列
185
186        Returns:
187            ndarray: （拟）价格序列
188        """
189        init_val = 1
190        # 从状态映射到涨跌幅。百分比。
191        price_map = {
192            0: -2,  # large drop
193            1: -1.5,  # drop
194            2: 0,  # no change
195            3: 1.5,  # small rise
196            4: 2  # large rise
197        }
198        enable_price_map = True
199        ret = np.zeros(len(output_list))
200        ret[0] = init_val
201        if enable_price_map:
202            for i in range(1, len(output_list)):
203                ret[i] = ret[i - 1] * (1 + price_map[output_list[i]] / 100)
204        else:
205            for i in range(1, len(output_list)):
206                ret[i] = ret[i - 1] + output_list[i] - 2
207        return ret
208    # 愚蠢的人类试图预测未来
209    n_seq = 28 # 预测长度
210    n_turn = 100 # 随机次数
211    
212    # == 使用前向算法，估计最后一个状态的概率分布 ==
213    # 由于已经在前面调用过，所以我们直接使用 alpha 计算结果
214    alpha_last = alpha[:, 99]
215    output_real = ob_seq[100:128]  # 实际值
216
217    # 最后状态。根据前向算法得出
218    forward_last_state = np.argmax(alpha_last)
219    price_seq_real = _reduce_output(output_real) # 真实价格序列
220    acc_seq = np.zeros(n_turn)
221    # 进行 n_turn 次模拟
222    for turn in range(n_turn):
223        last_state = forward_last_state
224        output_fit = np.zeros(n_seq)  # 预测值
225        for t in range(n_seq):
226            # 按照状态转移概率分布，抽取一个状态
227            state = np.random.choice(
228                np.arange(3), p=model.trans_pbb[last_state])
229            # 按照此状态到各个观测值的概率，抽取一个观测值
230            output = np.random.choice(np.arange(5), p=model.emit_pbb[state, :])
231            output_fit[t] = output
232
233            last_state = state
234
235        price_seq_fit = _reduce_output(output_fit)
236        # # 计算拟合优度
237        # # -- 残差平方和 RSS
238        # ss_res = np.sum((price_seq_real - price_seq_fit) ** 2)
239        # # -- 总平方和 TSS
240        # ss_tot = np.sum((price_seq_real - np.mean(price_seq_real)) ** 2)
241        # # -- 拟合优度，伪平方和
242        # r2 = 1 - (ss_res / ss_tot)
243        # acc_seq[turn] = r2
244        from seqeval.metrics import accuracy_score
245        acc_seq[turn] = accuracy_score(output_real, output_fit)
246        # acc_seq[turn] = accuracy_score(price_seq_real, price_seq_fit)
247        
248
249    # 输出拟合度的均值和标准方差
250    print("acc_seq=")
251    print(acc_seq)
252    print("acc_seq.mean=")
253    print(acc_seq.mean())
254    print("acc_seq.std=")
255    print(acc_seq.std())
256
257
258if __name__ == "__main__":
259    main()