目录
MCMC(一)蒙特卡罗方法
MCMC(二)马尔科夫链
MCMC(三)MCMC采样和M-H采样
MCMC(四)Gibbs采样
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 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]]
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为稳定概率分布。
下面的性质如何理解?
REF:
http://setosa.io/ev/ (可视化演示算法)