分类:Markov过程的阶和转移矩阵的计算
来自Big Physics
在博弈研究，以及其他很多研究中，Markov过程都是特别重要的概念。Markov过程就是指当前时刻的状态仅仅由有限步之前的状态决定，也就是，对于离散时间离散状态的系统来说，其当前时刻的决定机制可以表达成为一个矩阵，例如[math]\displaystyle{ P\left(S_{t}|S_{t-1}S_{t-2}...S_{t-L}\right) }[/math]。一阶Markov过程的意思就是[math]\displaystyle{ L=1 }[/math]。于是，如果系统有[math]\displaystyle{ N }[/math]个状态，则这个矩阵是[math]\displaystyle{ N\times N }[/math]的。更高阶的会表现为一个[math]\displaystyle{ N\times N^{L} }[/math]的矩阵。
在从实际数据中估计这个矩阵的过程中，有几个问题：
- 怎么知道哪一个[math]\displaystyle{ L }[/math]比较合适？
- 尤其在博弈研究中，牵涉到多个体参与者（这里的多个参与者仅仅是为了获得多个样本点，不是单次博弈所需要的人数。如果单次博弈实验需要人数是多人，则本来就要把这样的一群人当作系统来研究，系统的状态就是这群人合——其实是直积——起来的状态），每一个（或者每一群个体，这里的群就是单次实验本身所需要的参与者数量，例如2人博弈则需要2人）个体多轮的问题，估算的顺序可能会有影响，怎么决定：例如先通过每一个个体的时间序列来计算转移矩阵，然后对转移矩阵做群体平均；例如直接把所有人的时间序列拿来估算转移矩阵？前者称为按照过程样本来计算，后者称为按照系综来计算。
对于第一个问题，不能简单按照所估计出来的矩阵的理论解和实际分布函数的符合程度来决定。对于一个本质上就是两步的Markov过程，按照单步所估算出来矩阵，也能够给出和实际分布函数相符程度很高的理论分布函数（实际上，这个原因直观地说，就是在稳定状态，系统的高阶实际上就是相当于把其他各阶的信息用稳定分布折合回去。或者更极端地，可以写成一个“0-阶Markov过程”——仅仅用稳定分布来产生稳定分布。因此，仅仅看最终的分布函数，是不能当作检验的标准的）。当然，朴素来说，可以先按照高阶的计算出来，然后看看，是不是满足[math]\displaystyle{ P\left(S_{t}|S_{t-1}S_{t-2}\right)＝P\left(S_{t}|S_{t-1}\right) }[/math]来判断。更好的方法——效率更高、更定量的方法——应该其他人已经研究过了。需要做一些知识的补充。
对于第二个问题，按照分类:Surprised by the Gambler’s and Hot Hand Fallacies的结果，两种计算顺序有可能不一样。
下一步工作
- 看看Markov过程参数估计的方法（通常的条件概率确实是估计方法之一，在这里补充几篇文献），看看更好的确定[math]\displaystyle{ L }[/math]的方法，以及给出来这个估计出来的转移概率的符合程度（goodness of fit）的方法，例如[1][2][3][4]。尤其是[4] 。但是，也有一些方法显然是错的，例如[5]是不可靠的方法：高阶Markov的定义就有问题，参数估计方法也是有问题的。
- 看看在实际估算中，按照过程样本计算的，和，按照系综计算的，是否都有人采用。如果都有人采用，则，发出这样一个警告——两种计算可能给出来不同的转移矩阵，并且只有系综计算才是真的数学上的定义所对应的，是有学术贡献的。
参考文献
- ↑ Philipp Singer，Order estimation for Markov Chain Models, http://www.philippsinger.info/?p=44
- ↑ M.Papapetrou and D.Kugiumtzis, Markov chain order estimation with conditional mutual information, Physica A, 392(7), 2013, 1593-1601, https://doi.org/10.1016/j.physa.2012.12.017
- ↑ G. Morvai and B. Weiss Order estimation of Markov chains, IEEE Transactions on Information Theory 51(4), 1496 - 1497(2005)
- ↑ 4.0 4.1 Besag J, Mondal D. Exact goodness-of-fit tests for Markov chains. Biometrics. 69(2):488-96(2013). doi: 10.1111/biom.12009.
- ↑ Wai-Ki Ching, Michael K. Ng and Eric S. Fung, Higher-order multivariate Markov chains and their applications, Linear Algebra and its Applications 428 (2008) 492-507
附录
# http://www.bigphysics.org/index.php/%E5%88%86%E7%B1%BB:Markov%E8%BF%87%E7%A8%8B%E7%9A%84%E9%98%B6%E5%92%8C%E8%BD%AC%E7%A7%BB%E7%9F%A9%E9%98%B5%E7%9A%84%E8%AE%A1%E7%AE%97 #1-step Markovian tranfer matrix is defined as p1|0><0| + (1-p1)|1><0|+ (1-p2)|0><1| + p2|1><1| = [p1, 1-p2 \\1-p1, p2] #2-step Markovian tranfer matrix is defined as q00|0><00| + (1-q00)|1><00|+ q01|0><01| + (1-q01)|1><01| + q10|0><10| + (1-q10)|1><01|+ q11|0><11| + (1-q11)|1><11|= [q00, q01, q10, q11\\1-q00, 1-q01, 1-q10, 1-q11], |r_{t}><r_{t-1}r_{t-2}| import sys, getopt import random import numpy as np def GenerateSingleStep(r1, r2, M): x=random.uniform(0,1) #can be replaced with a better and specialized random number generator if x<M[r1*2+r2]: #M[r1,r2] r=0 else: r=1 return r def main(argv): L=100 #number of iterations N=10000 #length of sequences M=[0.5]*4 #value of elements in transfer matrix [M00, M01, M10, M11 \\ 1-M00, 1-M01, 1-M10, 1-M11] #take parameters from command-line input try: opts, args = getopt.getopt(argv, "h:L:N:", ["q00=","q01=","q10=","q11="]) except getopt.GetoptError: print ("Error: please use the command as Markov2.py -L <L> -N <N> --q00 <q00> --q01 <q01> --q10 <q10> --q11 <q11>") sys.exit(2) for opt, arg in opts: if opt == "-h": print("Markov2.py -L <L> -N <N> --q00 <q00> --q01 <q01> --q10 <q10> --q11 <q11>") sys.exit() elif opt =="-L": L = int(arg) elif opt == "-N": N = int(arg) elif opt =="--q00": M[0*2+0] = float(arg) elif opt == "--q01": M[0*2+1] = float(arg) elif opt =="--q10": M[1*2+0] = float(arg) elif opt == "--q11": M[1*2+1] = float(arg) #parameter input ends here P1=0;Q1=0; #final distribution P1ensemble = 0; Q1ensemble = 0 #number of 00 and 01 in the records of all time sequences P2ensemble = 0; Q2ensemble = 0 #number of 10 and 11 in the records of all time sequences records0=0; records1=0 #number of records of 0x and 1x in all time eqquences p1ensemble=0; p2ensemble=0 #final result of matrix over the whole ensemble p1sample=0; p2sample=0 #final result of matrix per trajectory first and then over all the trajectories for trial in range(L): P1sample=0; Q1sample=0; P2sample=0; Q2sample=0 #statistics per trajectory r1=random.randint(0, 1);r2=random.randint(0, 1) #value of r in the first two rounds for step in range(N-1): r=GenerateSingleStep(r1, r2, M) #use the transfer matrix P(r|r1,r2) to determine r for a given (r1,r2) #print(r2,r1,r) if r1==0: #in the case of r1==0 if r==0: P1sample=P1sample+1 #statistics per trajectory (trial) P1ensemble=P1ensemble+1 #statistics accumuated throughout the whole ensemble P1=P1+1 else: Q1sample=Q1sample+1 #statistics per trajectory (trial) Q1ensemble=Q1ensemble+1 #statistics accumuated throughout the whole ensemble Q1=Q1+1 else: #in the case of r1==1 if r==0: Q2sample=Q2sample+1 #statistics per trajectory (trial) Q2ensemble=Q2ensemble+1 #statistics accumuated throughout the whole ensemble P1=P1+1 else: P2sample=P2sample+1 #statistics per trajectory (trial) P2ensemble=P2ensemble+1 #statistics accumuated throughout the whole ensemble Q1=Q1+1 r2=r1 r1=r if P1sample>0 or Q1sample>0: #in the case of no records generated this round, one need this condition to avoid 0/0 records0=records0+1 p1sample=p1sample+1.0*P1sample/(1.0*P1sample+1.0*Q1sample) #matrix [p1, 1-p2 \\1-p1, p2] within a trajectory if P2sample>0 or Q2sample>0: #in the case of no records generated this round, one need this condition to avoid 0/0 records1=records1+1 p2sample=p2sample+1.0*P2sample/(1.0*P2sample+1.0*Q2sample) #matrix [p1, 1-p2 \\1-p1, p2] within a trajectory p1sample=p1sample/(1.0*records0); p2sample=p2sample/(1.0*records1) #matrix [p1, 1-p2 \\1-p1, p2] calculated for each sample sequences, then averaged over all the sequences p1ensemble=1.0*P1ensemble/(1.0*P1ensemble+1.0*Q1ensemble) p2ensemble=1.0*P2ensemble/(1.0*P2ensemble+1.0*Q2ensemble) #matrix [p1, 1-p2 \\1-p1, p2] directly calculated over all samples, thus from the whole ensemble Mensemble=np.array([[p1ensemble, 1.0-p2ensemble],[1.0-p1ensemble, p2ensemble]]) Eensemble,Pensemble=np.linalg.eig(Mensemble) Msample=np.array([[p1sample, 1.0-p2sample],[1.0-p1sample, p2sample]]) Esample,Psample=np.linalg.eig(Msample) print("L=", L, "N=", N, "q00=", M[0], "q01=", M[1], "q10=", M[2], "q11=", M[3]) print("transfer matrix calculated from ensemble:") print("[%5.3f, %5.3f \\\\ %5.3f, %5.3f]" %(Mensemble[0,0], Mensemble[0,1], Mensemble[1,0], Mensemble[1,1])) print("Stationary state of the theoretical transfer matrix is:") for i in range(2): if abs(Eensemble[i]-1.0)<1.0e-3: print("[%5.3f, %5.3f]" %(Pensemble[0,i]/(Pensemble[0,i]+Pensemble[1,i]), (Pensemble[1,i]/(Pensemble[0,i]+Pensemble[1,i])))) print("Real final distribution is:") print("[%5.3f, %5.3f]" %(1.0*P1/(1.0*P1+1.0*Q1), 1.0*Q1/(1.0*P1+1.0*Q1))) print("transfer matrix calculated from individual trajectories and then average all trjectories:") print("[%5.3f, %5.3f \\\\ %5.3f, %5.3f]" %(Msample[0,0], Msample[0,1], Msample[1,0], Msample[1,1])) print("Stationary state of the theoretical transfer matrix is:") for i in range(2): if abs(Esample[i]-1.0)<1.0e-3: print("[%5.3f, %5.3f]" %(Psample[0,i]/(Psample[0,i]+Psample[1,i]), (Psample[1,i]/(Psample[0,i]+Psample[1,i])))) if __name__ == "__main__": main(sys.argv[1:])
运行：
python3 Markov2.py -L 1 -N 100000 --q00 0.6 --q01 0.4 --q10 0.3 --q11 0.8
得到：
L= 1 N= 100000 q00= 0.6 q01= 0.4 q10= 0.3 q11= 0.8 transfer matrix calculated from ensemble: [0.501, 0.535 \\ 0.499, 0.465] Stationary state of the theoretical transfer matrix is: [0.517, 0.483] Real final distribution is: [0.517, 0.483] transfer matrix calculated from individual trajectories and then average all trjectories: [0.501, 0.535 \\ 0.499, 0.465] Stationary state of the theoretical transfer matrix is: [0.517, 0.483]
说明：这个实际上是一个二阶Markov过程，但是，按照统计出来的结果构造出来的一阶Markov过程的转移矩阵得到的稳定分布也和实际分布符合的非常好。甚至实际上，就可以变成0-阶Markov过程——也即是直接从分布函数中取样，也会符合得很好。因此，稳定分布是否符合不是确定阶数的好办法。
