分类: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]的矩阵。
在从实际数据中估计这个矩阵的过程中,有几个问题:
- 怎么知道哪一个[math]\displaystyle{ L }[/math]比较合适?
- 尤其在博弈研究中,牵涉到多个体参与者(这里的多个参与者仅仅是为了获得多个样本点,不是单次博弈所需要的人数。如果单次博弈实验需要人数是多人,则本来就要把这样的一群人当作系统来研究,系统的状态就是这群人合——其实是直积——起来的状态),每一个(或者每一群个体,这里的群就是单次实验本身所需要的参与者数量,例如2人博弈则需要2人)个体多轮的问题,估算的顺序可能会有影响,怎么决定:例如先通过每一个个体的时间序列来计算转移矩阵,然后对转移矩阵做群体平均;例如直接把所有人的时间序列拿来估算转移矩阵?前者称为按照过程样本来计算,后者称为按照系综来计算。
对于第一个问题,不能简单按照所估计出来的矩阵的理论解和实际分布函数的符合程度来决定。对于一个本质上就是两步的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的结果,两种计算顺序有可能不一样。
下一步工作
- 看看更好的确定[math]\displaystyle{ L }[/math]的方法
- 看看在实际估算中,按照过程样本计算的,和,按照系综计算的,是否都有人采用。
附录
# 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:])