参考：
【年终分享】彩票数据预测算法(一)：离散型马尔可夫链模型实现【附C#代码】…
例一
应用加权马尔可夫链模型预测某河流月份的降雨量
加权马尔可夫链模型预测的基础知识可以看这篇文章：
三种马尔可夫链预测方法
一、对历史数据进行分组
数据分组的方法：指标值序列分组方法
为了符合马尔可夫链预测模型的计算要求，把年降雨量分为5个等级。即5种状态，据调查分析，按照有关水文序列约定俗成的分法，
- 年降雨量小于450毫米为干旱年，定为状态1；
- 450~600毫米为偏旱年，定为状态2；
- 600~750毫米为正常年，定为状态3；
- 750~900毫米为偏丰年，定为状态4；
- 900毫米以上为丰涝年，定为状态5；。
以此标准划分各年降雨量的等级。如表（某河流5月份降雨量及其等级）所示：
一共22年的数据
二、“马氏性”检验
“马氏性”检验的知识：
马尔可夫性的统计检验（马氏性检验）
间隔1年，对这22年进行统计总结状态变化次数如下：
次数相加为n-1=22-1=21次（n是总的数据个数，也就是22年，1就是间隔为1年）
“边际概率”，用字母p·j表示 ，其中
此题的边际概率计算过程如下：
本次状态为1时，下次状态为1、2、3、4、5的频数为2+1+1+1+0=5次，而频数一共21次，那么状态为1的“边际概率”为P1j=5/21=0.238。同理可以计算状态2、3、4、5的“边际概率”为
当m很大时χ2（卡方分布） 统计量
（式子里面的log是因为Inx，在编程中Inx常写为logx）
如果Pij=0，Inx=-∞，所以对此处理是，当Pij=0，fij=0，那么0×∞=0。
它将服从自由度为(m−1)2的χ2（卡方）分布，现在给定显著性水平为α，经查表可得χ2α((m-1)2)的值 。如果χ2＞χ2α((m-1)2)，则拒绝零假设，可以认为序列具备“马氏性”，反之，则这个序列不能当作马尔可夫链来对待。
此题的卡方分布计算过程如下：
已知P1的矩阵为：
矩阵里的元素便是pij。
import numpy as np
m = 5 #状态数
f = np.array([[2,1,1,1,0], [2,3,3,0,0], [2,1,2,0,1], [0,1,0,0,0], [0,1,0,0,0]]) #频数矩阵
P = np.array([[2/5,1/5,1/5,1/5,0], [2/8,3/8,3/8,0,0], [2/6,1/6,2/6,0,1/6], [0,1,0,0,0], [0,1,0,0,0]]) #转移概率矩阵
Pj = np.array([5/21, 8/21, 6/21, 1/21, 1/21])
x = 0
for i in range(m):
for j in range(m):
if f[i][j] == 0:
x += 0
else:
x += f[i][j] * abs(np.log(P[i][j] / Pj[j]))
X_2 = 2*x
print('卡方分布X^2为：', X_2)
结果：
卡方分布X^2为： 18.850359569614124
给定显著性水平α=0.05，自由度（m-1）2=（5-1）2=16。查卡方分布的表格可知分位点X0.052=26.296。
由于本题计算的X2＜X0.052，故“年5月份降水量”所对应的的状态序列不满足“马氏性”，所以这个序列不能当作马尔可夫链来对待。
例二
改进的加权马尔可夫链预测方法及应用
加权马尔可夫链在降水状况预测中的应用
加权马尔可夫链预测模型的优化方法及应用
一、对历史数据进行分组
陕西省河曲水文站1952—1998年的年降水量数据见下表
import pandas as pd
rain = pd.read_excel(r'Rain.xlsx')
print(rain)
结果：
年份 降雨量（毫米）
0 1952 261.9
1 1953 486.4
2 1954 631.5
3 1955 259.0
4 1956 568.0
5 1957 398.2
6 1958 479.6
7 1959 697.6
8 1960 397.7
9 1961 640.4
10 1962 247.1
11 1963 387.7
12 1964 694.2
13 1965 211.4
14 1966 322.6
15 1967 656.6
16 1968 325.3
17 1969 603.8
18 1970 424.8
19 1971 383.3
20 1972 238.8
21 1973 423.0
22 1974 237.1
23 1975 330.7
24 1976 445.9
25 1977 518.9
26 1978 492.6
27 1979 490.3
28 1980 257.0
29 1981 400.0
30 1982 347.5
31 1983 363.8
32 1984 411.5
33 1985 356.2
34 1986 381.2
35 1987 317.0
36 1988 473.0
37 1989 373.7
38 1990 369.0
39 1991 348.3
40 1992 469.4
41 1993 228.1
42 1994 338.8
43 1995 546.1
44 1996 358.9
45 1997 412.0
46 1998 372.3
数据分组的方法：指标值序列分组方法
这里采用样本均值一均方差（标准差）分组法
对于数据序列x1，x2，……，xn，可看作是一个时间序列的前n个观测值，算出样本均值x和样本均方差s，根据具体情况以样本均值为中心，例如可将数据序列分成如下五组：
或：
import pandas as pd
import numpy as np
#以1952年——1995年年降水量为训练样本
rain = pd.read_excel(r'Rain.xlsx', skipfooter=3) #skipfooter:读取数据时，指定跳过的末尾行数。
print(rain)
N = rain.shape[0] #样本数量
print('样本数量：',N)
x_sum = sum(rain['降雨量（毫米）'])
x_Ave = x_sum/N #均值
print('均值x：',x_Ave)
s = np.sqrt(1/(N-1) * (sum((rain['降雨量（毫米）']-x_Ave)**2))) #均方差/标准差
print('均方差/标准差s：',s)
结果：
年份 降雨量（毫米）
0 1952 261.9
1 1953 486.4
2 1954 631.5
3 1955 259.0
4 1956 568.0
5 1957 398.2
6 1958 479.6
7 1959 697.6
8 1960 397.7
9 1961 640.4
10 1962 247.1
11 1963 387.7
12 1964 694.2
13 1965 211.4
14 1966 322.6
15 1967 656.6
16 1968 325.3
17 1969 603.8
18 1970 424.8
19 1971 383.3
20 1972 238.8
21 1973 423.0
22 1974 237.1
23 1975 330.7
24 1976 445.9
25 1977 518.9
26 1978 492.6
27 1979 490.3
28 1980 257.0
29 1981 400.0
30 1982 347.5
31 1983 363.8
32 1984 411.5
33 1985 356.2
34 1986 381.2
35 1987 317.0
36 1988 473.0
37 1989 373.7
38 1990 369.0
39 1991 348.3
40 1992 469.4
41 1993 228.1
42 1994 338.8
43 1995 546.1
样本数量： 44
均值x： 414.4318181818181
均方差/标准差s： 130.3853443981594
上面的错误改为：
这里α1与α4在划分区域对称，应该是相同的；同理，α2与α3相等。
import pandas as pd
import numpy as np
#以1952年——1995年年降水量为训练样本
rain = pd.read_excel(r'Rain.xlsx', skipfooter=3) #skipfooter:读取数据时，指定跳过的末尾行数。
N = rain.shape[0] #样本数量
x_sum = sum(rain['降雨量（毫米）'])
x_Ave = x_sum/N #均值
s = np.sqrt(1/(N-1) * (sum((rain['降雨量（毫米）']-x_Ave)**2))) #均方差/标准差
#分级标准
#(-∞,x-a1*s)，(x-a1*s,x-a2*s)，(x-a2*s,x+a3*s)，(x+a3*s,x+a4*s)，(x+a4*s,+∞) #x为样本均值，s为样本均方差/标准差
a1, a2, a3, a4 = 1.1, 0.3, 0.3, 1.1
#状态1为丰，状态2为偏丰，状态3为正常，状态4为偏枯，状态5为枯
rain_class_list = []
for i in range(N):
if rain['降雨量（毫米）'][i] >= x_Ave+a4*s :
rain_class_list.append(1)
elif rain['降雨量（毫米）'][i] >= x_Ave+a3*s :
rain_class_list.append(2)
elif rain['降雨量（毫米）'][i] >= x_Ave-a2*s :
rain_class_list.append(3)
elif rain['降雨量（毫米）'][i] >= x_Ave-a1*s :
rain_class_list.append(4)
else:
rain_class_list.append(5)
print(x_Ave+a1*s, x_Ave+a2*s, x_Ave-a2*s, x_Ave-a1*s)
# rain['状态'] = pd.Series(rain_class_list)
rain['状态'] = rain_class_list
print(rain)
结果：
557.8556970197934 453.5474215012659 375.31621486237026 271.00793934384274
年份 降雨量（毫米） 状态
0 1952 261.9 5
1 1953 486.4 2
2 1954 631.5 1
3 1955 259.0 5
4 1956 568.0 1
5 1957 398.2 3
6 1958 479.6 2
7 1959 697.6 1
8 1960 397.7 3
9 1961 640.4 1
10 1962 247.1 5
11 1963 387.7 3
12 1964 694.2 1
13 1965 211.4 5
14 1966 322.6 4
15 1967 656.6 1
16 1968 325.3 4
17 1969 
最低0.47元/天 解锁文章
1277
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
