概率机器人 第二章第2道习题 马尔可夫链 答案

#coding=utf-8

# 用于按概率生成随机数
import numpy as np

# 函数；生成未来若干天的天气
# now : 今天天气， 0-晴天  1-多云  2-雨
# count： 返回多少天之后的天气的序列，默认只返回第二天
# trans: 转移矩阵
# 返回：天气list
def gen_weather( now, count, trans):
	weather_list = []
	for num in range( count ):
		# 根据今天天气对应的明天天气的转换概率，生成第二天的天气
		p = np.array( trans[ now] )
		#   0-晴天  1-多云 2-雨
		now = np.random.choice([0,1,2], p=p.ravel())
		# 保存到天气链表
		weather_list.append( now )
	return weather_list

# 函数：计算马尔可夫链的平稳分布(numpy array 方式）
# init_array: 平稳分布的初始值list, 按[晴，多云，雨]顺序
# transfer_mat: 转移矩阵
# error: 收敛的误差值阈值，如1e-07
# 返回： 平稳分布list
def np_markov(init_array, transfer_mat, error):
	restmp = init_array
	count = 0
	while True:
		#print ( count, restmp )
		res = np.dot( restmp, transfer_mat )
		diff = res - restmp
		if diff.max() < error:
			break;
		restmp = res
		count += 1
	return count, restmp

# 函数： 计算马尔可夫链的平稳分布(python array 方式）
# init_array: 平稳分布的初始值list, 按[晴，多云，雨]顺序
# trans: 转移矩阵
# error: 收敛的误差值阈值，如1e-07
# 返回： 平稳分布list
def array_markov( init_array, trans, error ):
	count = 0  # 计数器
	while True:
		#print ( count, init_array )
		# 计算 w = init_array*trans
		w = [0,0,0]
		# 计算init_array*trans中第二天为晴的概率
		w[0] = init_array[0] * trans[0][0] + init_array[1]*trans[1][0] + init_array[2]*trans[2][0]
		# 计算init_array*trans中第二天为多云 的概率
		w[1] = init_array[0] * trans[0][1] + init_array[1]*trans[1][1] + init_array[2]*trans[2][1]
		# 计算init_array*trans中第二天为雨 的概率
		w[2] = init_array[0] * trans[0][2] + init_array[1]*trans[1][2] + init_array[2]*trans[2][2]

		# 得到init_array与w的delta
		delta = [0,0,0]
		for i in range(3):
			delta[i] = abs( init_array[i] - w[i])
		# 刷新init_array, 进行下一轮计算
		init_array = w
		# 若差值最够小，则停止计算
		if delta[0] < error and delta[1] < error and delta[2] < error :
			break;
		count += 1
	return count, w


# 主程序
if __name__ == "__main__":
	# 天气状态转换矩阵
	#   第一列/行：晴
	#   第二列/行：多云
	#   第三列/行：雨
	trans = [
	   [0.8, 0.2, 0.0],
	   [0.4, 0.4, 0.2],
	   [0.2, 0.6, 0.2]
	]

	# 开始生成随机天气
	print ("\n根据状态转换矩阵生成随机天气")
	start_weather = 1   #起始天气，0-晴天  1-多云  2-雨
	count = 20    # 生成多少天的天气
	print ( "起始天气:%d, 生成天数:%d "% (start_weather,count) )
	li  = gen_weather( start_weather, count, trans )
	print ( "计算得到天气序列为: ",li )

	# 计算马尔可夫链的平稳分布(方式1：使用numpy)
	# 天气状态转换矩阵
	#   第一列/行：晴
	#   第二列/行：多云
	#   第三列/行：雨
	transfer_mat = np.array(
	  [[0.8, 0.2, 0.0],
	   [0.4, 0.4, 0.2],
	   [0.2, 0.6, 0.2]]
	)
	# np的精度是小数点后16位
	np.set_printoptions(precision=16)
	print ("\n计算马尔可夫链的平稳分布(方式1):")
	# 设置初始分布init_array(任意){晴，多云，雨}
	init_array = np.array([0.3, 0.4, 0.3 ])
	# 停止计算的误差值
	error = 1e-14
	print ( "初始分布:",init_array, " error:", error)
	count, li = np_markov(init_array, transfer_mat, error )
	print ("计算次数:", count, "\n平稳分布为: ", li )



	# 计算马尔可夫链的平稳分布(方式2：使用python array)
	print ("\n计算马尔可夫链的平稳分布(方式2):")
	# 设置初始分布init_array(任意){晴，多云，雨}
	init_array = [0.3, 0.4, 0.3 ]
	# 停止计算的误差值
	error = 1e-14
	print ( "初始分布:",init_array, " error:", error)
	count, li = array_markov( init_array ,trans, error )
	print ("计算次数:", count, "\n平稳分布为: ", li )

运行结果：

根据状态转换矩阵生成随机天气
起始天气:1, 生成天数:20
计算得到天气序列为:  [0, 0, 0, 0, 0, 0, 1, 2, 1, 1, 1, 2, 0, 1, 0, 1, 0, 0, 1, 2]

计算马尔可夫链的平稳分布(方式1):
初始分布: [0.3 0.4 0.3]  error: 1e-14
计算次数: 43
平稳分布为:  [0.6428571428571346 0.2857142857142918 0.0714285714285755]

计算马尔可夫链的平稳分布(方式2):
初始分布: [0.3, 0.4, 0.3]  error: 1e-14
计算次数: 43
平稳分布为:  [0.6428571428571395, 0.2857142857142889, 0.07142857142857345]