分类:Markov过程的阶和转移矩阵的计算

在博弈研究，以及其他很多研究中，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]的矩阵。

在从实际数据中估计这个矩阵的过程中，有几个问题：

1. 怎么知道哪一个[math]\displaystyle{ L }[/math]比较合适？
2. 尤其在博弈研究中，牵涉到多个体参与者（这里的多个参与者仅仅是为了获得多个样本点，不是单次博弈所需要的人数。如果单次博弈实验需要人数是多人，则本来就要把这样的一群人当作系统来研究，系统的状态就是这群人合——其实是直积——起来的状态），每一个（或者每一群个体，这里的群就是单次实验本身所需要的参与者数量，例如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的结果，两种计算顺序有可能不一样。

下一步工作

1. 看看更好的确定[math]\displaystyle{ L }[/math]的方法，例如^[1]
2. 看看在实际估算中，按照过程样本计算的，和，按照系综计算的，是否都有人采用。如果都有人采用，则，发出这样一个警告——两种计算可能给出来不同的转移矩阵，并且只有系综计算才是真的数学上的定义所对应的，是有学术贡献的。

参考文献

引用错误：在<references>中以“Papapetrou”名字定义的<ref>标签没有在先前的文字中使用。
引用错误：在<references>中以“Morvai”名字定义的<ref>标签没有在先前的文字中使用。

附录

# 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; P2ensemble = 0; Q2ensemble = 0 #number of 1s in the records of all time sequences
	p1ensemble=0; p2ensemble=0 #final result of matrix per ensemble
	for trial in range(L):
		P1sample=0; Q1sample=0; P2sample=0; Q2sample=0 #statistics per sequence
		r1=0;r2=0 #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) 
			if r1==0:                       #in the case of r1==0
				P1=P1+1
				if r==0:
					P1ensemble=P1ensemble+1 #statistics accumuated throughout the whole ensemble
				else:
					Q1ensemble=Q1ensemble+1 #statistics accumuated throughout the whole ensemble
			else:                           #in the case of r1==1
				Q1=Q1+1
				if r==0:
					P2ensemble=P2ensemble+1 #statistics accumuated throughout the whole ensemble
				else:
					Q2ensemble=Q2ensemble+1 #statistics accumuated throughout the whole ensemble
			r2=r1
			r1=r

	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
	
	M2=np.array([[p1ensemble, p2ensemble],[1.0-p1ensemble, 1.0-p2ensemble]])
	E,P=np.linalg.eig(M2)
	
	print("L=", L, "N=", N, "q00=", M[0], "q01=", M[1], "q10=", M[2], "q11=", M[3])
	print("tranfer matrix calculated from ensemble:")
	print("[%5.3f, %5.3f \\\\ %5.3f, %5.3f]" %(M2[0,0], M2[0,1], M2[1,0], M2[1,1]))
	print("Stationary state of the theoretical transfer matrix is:")
	for i in range(2):
		if abs(E[i]-1.0)<1.0e-3:
			print("[%5.3f, %5.3f]" %(P[0,i]/(P[0,i]+P[1,i]), (P[1,i]/(P[0,i]+P[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)))
	
	

if __name__ == "__main__":
    main(sys.argv[1:])