马尔科夫链预测，Python实现

–记录自己学习的步伐，点滴的生活，以后学习和复习使用。纯手打，代码不抄袭。–

来源百度百科，具体定义和性质可以查看百度百科的内容。马尔科夫链
或者看一下下面的这个文章，也是翻译转载过来的。偷渡一下：[译] 用 Python 实现马尔可夫链的初级教程

用概率数学公式表示如下：
Pr( Xn+1 = x | X1 = x1, X2 = x2, …, Xn = xn) = Pr( Xn+1 = x | Xn = xn)

就是说 Xn+1 的概率只和之前的 Xn 的概率有关。所以只需要知道上一个状态就可以确定现在状态的概率分布，满足条件独立（也就是说：只需要知道现在状态就可以确定下一个状态，而之前的所有状态，都不会影响）。

例1：
某计算机机房的一台计算机经常出现故障，研究者每隔15分钟观察一次计算机的运行状态，手机了24小时的数据，用1表示正常，0表示不正常。数据如下：

str1 =
"1110010011111110011110111111001111111110001101101111011011010111101110111101111110011011111100111"

假定这就是无规律，白噪声的序列，我也没有做过白噪声检验，反正就用这个例子做了。
现在有四种情况，分别是：

第一步1之后出现下一步0，表示为10
第一步1之后出现1，表示为11
0之后出现0，表示为00
0之后出现1，表示为01

那么要预测下一个时间点将会出现的状态。
首先要计算上面4中情况一共出现多少次，然后才能计算概率。
比如111，里面出现11的情况有两次，而不是一次。这就要做循环统计了。
统计出每一个多少，然后做转移矩阵。样式跟下面这个类似

	0	1
0	…	…
1	…	…

其中index上面的0表示当前的一步，就第一步，columns上面的0表示下一步。也不是绝对的要这种形式，只要自己看得懂就行，自己有个习惯，并记住也能让别人看懂就好。


# 计算ab在s中出现的次数
def find_count(s,ab):ab_sum = 0for i in range(0,len(s)-1):if s[i:i+2] == ab:ab_sum+=1return ab_sum# 转移矩阵
def str_count_df(s):# 获得里面不重复的元素unique_items = np.unique(list(s))# 获得不重复元素个数n = unique_items.size# 默认行是这一次的，列是下一次的。类容是他们的转换情况df_ = pd.DataFrame(index=unique_items,columns=unique_items)for i in unique_items:for j in unique_items:df_.loc[i,j] = find_count(s,i+j)return df_# 转移矩阵，概率
def str_count_df_p(s):# 获得里面不重复的元素unique_items = np.unique(list(s))# 获得不重复元素个数n = unique_items.size# 默认行是这一次的，列是下一次的。类容是他们的转换情况df_ = pd.DataFrame(index=unique_items,columns=unique_items)for i in unique_items:for j in unique_items:df_.loc[i,j] = find_count(s,i+j)df_ = df_.div(df_.sum(axis=1),axis='index')return df_

上面两个转移矩阵代码是一样的，只不过最后返回值不一样，一个是转换成概率了，一个没转。
用上面的函数，将str1带进去就会得到转移矩阵

表示0→1共发生18次，以此类推…

表示如果当前是0，那么下一个是0的概率为30.77%，下一步为1的概率为69.23%。
每一行的和是1.这好理解的。
对当前数据，最后一个为1，那么预测下一步，有74.3%的概率不发生故障。

上面的方法不仅仅局限于两个类别的0和1，任何类别的转换都可以。
比如下面的观测系统数据：

str2 = "4321431123212344331113321222442323112431"

带入也可直接得到转移矩阵

这计算的都是一步转移概率矩阵。

如果是两步转移矩阵，就是一步转移矩阵的二次方。三步转移矩阵就是三次方，n步就是n次方

这是数学上已经证明过的。

现在我们假设下一步出现1,2,3,4的概率为 [0.1,0.3,0.2,0.4]
那么第二步出现的概率就可以算出来，用第一次的乘以转移矩阵：

#上面的转移矩阵为 df_2
np.array([0.1,0.3,0.2,0.4]).dot(np.array(df_2,dtype=np.float64))

out:

array([0.19454545, 0.22441558, 0.38402597, 0.19701299])

那么第三步：就是第一步*df_2的二次方,以此类推。
当n趋于无穷，最终的结果会不变。

# 这是第四步
np.array([0.1,0.3,0.2,0.4]).dot(np.linalg.matrix_power(np.array(df_2,dtype=np.float64),3))

out:

array([0.29265661, 0.28825609, 0.27010259, 0.14898471])

# 这是第10步
np.array([0.1,0.3,0.2,0.4]).dot(np.linalg.matrix_power(np.array(df_2,dtype=np.float64),9))

out:

array([0.29393371, 0.28893058, 0.26829268, 0.14884303])

# 这是第20步
np.array([0.1,0.3,0.2,0.4]).dot(np.linalg.matrix_power(np.array(df_2,dtype=np.float64),19))

out:

array([0.29393371, 0.28893058, 0.26829268, 0.14884303])

通过上面可以看出，最后已经不再变化。
不管你一开始的概率，也就是第一步的概率是多少，最后都会是这个结局。
跟第一步没关系。
这就是这个转移矩阵的极限分布。
这个在市场分析，或者市场占有率上，很有用。最后趋于一个恒定值，就是拥有多大的市场份额。
（假如这1,2,3,4分别代表不同的商品的话）

另外这个极限值是可以求出来的。我们假设它为（p1,p2,p3,p4）
可以知道这个向量，乘以转移矩阵以后还是它这个向量本身，所以不管它乘以多少个转移矩阵，都不会变化。
假如转移矩阵为:

[[a1,a2,a3,a4]
[b1,b2,b3,b4]
[c1,c2,c3,c4]
[d1,d2,d3,d4]]

根据矩阵乘积运算。有：

p1 X a1+p2Xb1+p3Xc1+p4Xd1=p1
p1 X a2+p2Xb2+p3Xc2+p4Xd2=p2
p1 X a3+p2Xb3+p3Xc3+p4Xd3=p3
p1 X a4+p2Xb4+p3Xc4+p4Xd4=p4

就是1 * 4的矩阵乘以4 * 4的矩阵。等于1 * 4的矩阵。

并且p1+p2+p3+p4=1

通过这5个方程就可以求出来p1,p2,p3,p4

用Python计算的时候，如果用scipy.linalg里面的solve，会有问题，因为这个方法要求传入的是一个方阵，可是这里有五个方程，四个未知数，无法求出。
我试过很多次，都不行，但可以去掉一个方程求得近似解。

我也是找了很久，有这么一个方法可以求：

# Return the least-squares solution to a linear matrix equation.
# 用最小二乘法来求解，a是系数矩阵，b是结果
np.linalg.lstsq(a,b)

代码：

# 注意看上面的方程，原转移矩阵要转置，并减去一个单位阵。
a = np.array(df_2,dtype=np.float64).T
a = a-np.eye(4)
a = np.append(a,np.array([1,1,1,1])).reshape(-1,4)
np.linalg.lstsq(a,np.array([0,0,0,0,1]))

out:

# 第一个返回值为我们要求得p1到p4
(array([0.29393371, 0.28893058, 0.26829268, 0.14884303]),array([1.02992113e-33]),4,array([2.01437748, 1.25270383, 1.05972215, 0.70596599]))
# 第二个返回值为残差。第三个返回值时秩，第四个为矩阵奇异值

并不是每一个转移矩阵n次方后都会平稳趋于一个定矩阵，这个平稳等性质，可以搜别的文章，有的不趋于，有的单数次趋于一个偶数次趋于一次。等等。

如果对一个转移概率矩阵进行n次方，他得到的值越来越小，但是每个值占的百分比却一直不变，那有可能是数值的精度问题，碰到过，后来转成分数就好了。或者求一下极限分布。

下面看一个有意思的问题，我也是在知乎上看到的。
不过他只给出了四个格子的步骤和方法。

如下图所示的迷宫共有9个格子，相邻格子有门相通，9号格子就是迷宫出口.
整个迷宫将会在5分钟后坍塌. 1号格子有一只老鼠，这只老鼠以每分钟一格的速度在迷宫里乱窜（它通过各扇门的机会均等）。
求此老鼠在迷宫坍塌之前逃生的概率。如果这只老鼠速度提高一倍，则老鼠在迷宫坍塌之前逃生的概率能增加多少？

以上是2018年上海应用数学知识竞赛的一道初赛试题。

对于这个问题，先求出他的转移概率矩阵，那么一切就好办了。
9 * 9的矩阵，而且从1到2和3的概率相同，其他也一样。如下：

# 写出每个格子的转移矩阵
a = np.array([[0,0.5,0,0.5,0,0,0,0,0],[1/3,0,1/3,0,1/3,0,0,0,0],[0,0.5,0,0,0,0.5,0,0,0],[1/3,0,0,0,1/3,0,1/3,0,0],[0,0.25,0,0.25,0,0.25,0,0.25,0],[0,0,1/3,0,1/3,0,0,0,1/3],[0,0,0,0.5,0,0,0,0.5,0],[0,0,0,0,1/3,0,1/3,0,1/3],[0,0,0,0,0,0,0,0,1]
])

由于第一步在1的位置，所以第一步走下去的概率为：

f = np.array([0,0.5,0,0.5,0,0,0,0,0])

那么计算5步的时候概率：

f.dot(np.linalg.matrix_power(a,4))

out:
array([0. , 0.27777778, 0. , 0.27777778, 0. ,0.16666667, 0., 0.16666667, 0.11111111])
所以停在9的概率也就是出去的概率为0.111111111

速度加快就是一分钟走两步，五分钟10步。

f.dot(np.linalg.matrix_power(a,9))

array([0.1399177 , 0. , 0.11316872, 0. , 0.22633745,
0. , 0.11316872, 0. , 0.40740741])

假定时间趋于很大的时候，求它的极限分布：
（代码不重复写了。根据那个公式构造）

array([ 1.42076747e-17, -5.55111512e-17, -1.66533454e-16, -8.67361738e-17,-3.33066907e-16, -1.24900090e-16, -1.11022302e-16, -6.93889390e-17,1.00000000e+00])

可以看得出来，最后9的概率等于1，其他的等于0.所以最终停留的地方一定是9，也就是出口。

继续，
倘若迷宫不存在出口，并且每个门通过概率一样。那么他最终停在每个格子上的概率是不是相同的？
重新构造他的转移概率矩阵：（其实就是最后一个行有变化，其他都一样）

a = np.array([[0,0.5,0,0.5,0,0,0,0,0],[1/3,0,1/3,0,1/3,0,0,0,0],[0,0.5,0,0,0,0.5,0,0,0],[1/3,0,0,0,1/3,0,1/3,0,0],[0,0.25,0,0.25,0,0.25,0,0.25,0],[0,0,1/3,0,1/3,0,0,0,1/3],[0,0,0,0.5,0,0,0,0.5,0],[0,0,0,0,1/3,0,1/3,0,1/3],[0,0,0,0,0,0.5,0,0.5,0]
])

求第103步和104步的概率分布情况

f = np.array([0,0.5,0,0.5,0,0,0,0,0])
print(f.dot(np.linalg.matrix_power(a,102)))
print(f.dot(np.linalg.matrix_power(a,103)))

out：

[0，0.25，0.，0.25，0， 0.25，0，0.25，0]
[0.16666667，0.        ， 0.16666667 ，0.        ， 0.33333333 0.0.16666667 ，0.         ，0.16666667]

多算几个你就会发现他是上面两组的交替，也就是奇数、偶数次会有不同的结果。

我们用上面公式求得他的极限分布p1-p9：

a_ = a.T-np.eye(9)
a_ = np.append(a_,np.array([1,1,1,1,1,1,1,1,1])).reshape(-1,9)
np.linalg.lstsq(a_,np.array([0,0,0,0,0,0,0,0,0,1]))

array([0.08333333, 0.125     , 0.08333333, 0.125     , 0.16666667,0.125     , 0.08333333, 0.125     , 0.08333333])

结果其实就是上面交替的两个矩阵的平均值。这就是他们最终每个点停靠的概率。

马尔科夫链预测，Python实现相关推荐

马尔可夫预测matlab编程,MATLAB马尔科夫链预测法
MATLAB马尔科夫链预测法发布时间:2018-06-10 09:10, 浏览次数:1195 , 标签: MATLAB 1,确定初始概率分布p0 输入初始状态矩阵: >>A=[11 19 ...
灰色马尔科夫链matlab,基于灰色-马尔科夫模型的电力功率预测
利用1998-2009每年的用电量预测2010年的用电量 QQ图片20130515210109.jpg (20.32 KB, 下载次数: 18) 1998-2009每年用电量数据 2013-5-15 ...
基于python的马尔科夫链在股价预测中的应用（基于Tushare）
TushareID:503535 文章目录前言一.马尔科夫链是什么二.代码如下三.马尔科夫预测模型在股价预测中的应用 1.数据来源及状态划分 2.状态转移概率矩阵 P 及初始状态概率向量图片 ...
通过Python实现马尔科夫链蒙特卡罗方法的入门级应用
通过把马尔科夫链蒙特卡罗(MCMC)应用于一个具体问题,本文介绍了 Python 中 MCMC 的入门级应用. GitHub 地址:https://github.com/WillKoehrsen/ai ...
python做马尔科夫模型预测法_通过Python的Networkx和Sklearn来介绍隐性马尔科夫模型...
Python部落(python.freelycode.com)组织翻译,禁止转载,欢迎转发. 文章梗概马尔科夫是何人? 马尔科夫性质是什么? 马尔科夫模型是什么? 是什么让马尔科夫模型成为隐性的? ...
【深度剖析HMM（附Python代码）】1.前言及隐马尔科夫链HMM的背景
1. 前言隐马尔科夫HMM模型是一类重要的机器学习方法,其主要用于序列数据的分析,广泛应用于语音识别.文本翻译.序列预测.中文分词等多个领域.虽然近年来,由于RNN等深度学习方法的发展,HMM模型逐 ...
用Python中的马尔科夫链进行营销渠道归因
用Python中的马尔科夫链进行营销渠道归因 --第一部分:"更简单 "的方法任何积极开展营销活动的企业都应该对确定哪些营销渠道推动了实际转化率感兴趣.投资回报率(ROI)是一个 ...
预测未来的神技——有趣的马尔科夫链
本文为量子金服原创文章,转载须授权作为概率论的一个重要分支,随机过程撑起了概率论的半壁江山,如今,它广泛使用在诸如天气预报.统计物理.天体物理.运筹决策.经济数学.安全科学.人口理论.可靠性及计算机 ...
python：马尔科夫链
参考学习网址:MCMC(二)马尔科夫链 - 刘建平Pinard - 博客园 (cnblogs.com) 目录一.什么是马尔科夫链 1.介绍 2.基本概念 3.马尔科夫链模型分析二.马尔科夫链模型状 ...

马尔科夫链预测，Python实现

马尔科夫链预测，Python实现相关推荐

最新文章

热门文章