目录

MCMC(一)蒙特卡罗方法

MCMC(二)马尔科夫链

MCMC(三)MCMC采样和M-H采样
MCMC(四)Gibbs采样 

 

 

马尔科夫链_Python

 

马尔科夫链_html_02

 

马尔科夫链_f5_03

 

马尔科夫链_Python_04

 

马尔科夫链_Python_05

 Python 2.7 版本:

import numpy as np
matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float)
vector1 = np.matrix([[0.3,0.4,0.3]], dtype=float)
for i in range(100):
    vector1 = vector1*matrix
    print "Current round:" , i+1
    print vector1

Python 3.7 版本:

import numpy as np
matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float) 
vector1 = np.matrix([[0.3,0.4,0.3]], dtype=float)
for i in range(100):
    vector1 = vector1*matrix
    print("Current round:" , i+1)
    print(vector1)

部分输出结果如下:

Current round: 1
[[ 0.405   0.4175  0.1775]]
Current round: 2
[[ 0.4715   0.40875  0.11975]]
Current round: 3
[[ 0.5156  0.3923  0.0921]]
Current round: 4
[[ 0.54591   0.375535  0.078555]]
。。。。。。
Current round: 58
[[ 0.62499999  0.31250001  0.0625    ]]
Current round: 59
[[ 0.62499999  0.3125      0.0625    ]]
Current round: 60
[[ 0.625   0.3125  0.0625]]
。。。。。。
Current round: 99
[[ 0.625   0.3125  0.0625]]
Current round: 100
[[ 0.625   0.3125  0.0625]]

马尔科夫链_Python_06

  Python 2.7 版本:

matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float)
vector1 = np.matrix([[0.7,0.1,0.2]], dtype=float)
for i in range(100):
    vector1 = vector1*matrix
    print "Current round:" , i+1
    print vector1

Python 3.7 版本:

import numpy as np
matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float)
vector1 = np.matrix([[0.7,0.1,0.2]], dtype=float)
for i in range(100):
    vector1 = vector1*matrix
    print("Current round:" , i+1)
    print(vector1)

 部分输出结果如下:

Current round: 1
[[ 0.695   0.1825  0.1225]]
Current round: 2
[[ 0.6835   0.22875  0.08775]]
Current round: 3
[[ 0.6714  0.2562  0.0724]]
Current round: 4
[[ 0.66079   0.273415  0.065795]]
。。。。。。。
Current round: 55
[[ 0.62500001  0.31249999  0.0625    ]]
Current round: 56
[[ 0.62500001  0.31249999  0.0625    ]]
Current round: 57
[[ 0.625   0.3125  0.0625]]
。。。。。。。
Current round: 99
[[ 0.625   0.3125  0.0625]]
Current round: 100
[[ 0.625   0.3125  0.0625]]

马尔科夫链_f5_07

  Python 2.7 版本:

matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float)
for i in range(10):
    matrix = matrix*matrix
    print "Current round:" , i+1
    print matrix

Python 3.7 版本:

import numpy as np
matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float)
for i in range(10):
    matrix = matrix*matrix
    print("Current round:" , i+1)
    print(matrix)

 输出结果如下:

Current round: 1
[[ 0.8275   0.13375  0.03875]
 [ 0.2675   0.66375  0.06875]
 [ 0.3875   0.34375  0.26875]]
Current round: 2
[[ 0.73555   0.212775  0.051675]
 [ 0.42555   0.499975  0.074475]
 [ 0.51675   0.372375  0.110875]]
。。。。。。
Current round: 5
[[ 0.62502532  0.31247685  0.06249783]
 [ 0.6249537   0.31254233  0.06250397]
 [ 0.62497828  0.31251986  0.06250186]]
Current round: 6
[[ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]]
Current round: 7
[[ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]]
。。。。。。
Current round: 9
[[ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]]
Current round: 10
[[ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]
 [ 0.625   0.3125  0.0625]]  

X*Y...Y=Z  ==> Y...Y=Z, Y和Z确定之后,X为任意矩阵都可以。Y的n次幂=Z 
Y为状态转移矩阵;Z为稳定概率分布。

 

马尔科夫链_html_08

下面的性质如何理解?

 

马尔科夫链_Python_09

 

马尔科夫链_html_10

 

马尔科夫链_f5_11

 

马尔科夫链_html_12

 


REF:


https://mp.weixin.qq.com/s?__biz=MzU1NTUxNTM0Mg==&mid=2247489378&idx=1&sn=4507090d90f09561e5a675a6d4106893&chksm=fbd27bc3cca5f2d581114a5f2f587d58fc57620399d9b909b8c07a1d152f685b782aab03fbf5&mpshare=1&scene=1&srcid=11273P7xL2PHQAcb1mgI9BDj#rd

http://setosa.io/ev/  (可视化演示算法)