python HiKat1创建的sdm0711-with-graph.py - https://repl.it/J59N/1

Posted 2021-05-09

#!/usr/bin/env python
# coding:utf-8
#
# use Python3 !!!
#
import math
from scipy.stats import norm

b1 = 1
b2 = 0
b3 = 2
b4 = 8

x = [0]
x.append(10 * b1 + 30)
x.append(10 * b2 + 30)
x.append(10 * b3 + 30)
x.append(10 * b4 + 30)
print("Your x")
print(x)


sigma = 10
gamma = 20
print("sigma:")
print(sigma)
print("gamma:")
print(gamma)

Mu = [50] 
P = [100]
Kalman = [0]
# P1|0 = P0, P2|1 = l
P_t_tmin1 = [-99999, 100]

# P_t|t-1を求める関数
def calc_P_t_tmin1(t):
	return P[t-1] + gamma

def calc_KalmanGain_at_t(t):
	return (P_t_tmin1[t]) / (P_t_tmin1[t] + sigma)

def calc_P_at_t(t):
	return (1.0 - Kalman[t]) * (P_t_tmin1[t])

def calc_Mu_at_t(t):
	return Mu[t-1] + Kalman[t] * (x[t] - Mu[t-1])


def calc_K_P_Mu_iter():
	for t in range(1,5):
		if t != 1:
			# Pt|t-1を計算して追加
			P_t_tmin1.append(calc_P_t_tmin1(t))

		# Kalman Gainを計算して追加
		Kalman.append(calc_KalmanGain_at_t(t))

		# Pを計算して追加
		P.append(calc_P_at_t(t))

		# Muを計算して追加
		Mu.append(calc_Mu_at_t(t))


calc_K_P_Mu_iter()
print("Mu:")
print(Mu)
print("P")
print(P)
print("Kalman Gain")
print(Kalman)
print("P_t|t-1")
print(P_t_tmin1)

P_x1 = norm.pdf(x = x[1] ,loc = Mu[0], scale = math.sqrt(P_t_tmin1[1] + sigma))
ln_P_x1 = math.log(P_x1)
print("\n")
print("Answer Q1 ================================")
print("lnP(x1)=")
print(ln_P_x1)
print("P(z1|x1)=")
print("N(z1;" + str(Mu[1]) + "," + str(P[1]) + ")")
print("\n")


P_x2_given_x1 = norm.pdf(x = x[2] ,loc = Mu[1], scale = math.sqrt(P_t_tmin1[2] + sigma))
ln_P_x2_given_x1 = math.log(P_x2_given_x1)
print("Answer Q2 ================================")
print("P(z2|x1)=")
print("N(z2;" + str(Mu[1]) + "," + str(P_t_tmin1[2]) + ")")
print("lnP(x2|x1)=")
print(ln_P_x2_given_x1)
print("lnP(x1,x2)=")
ln_P_x1_x2 = ln_P_x2_given_x1 + ln_P_x1
print(ln_P_x1_x2)
print("P(z2|x1,x2)=")
print("N(z2;" + str(Mu[2]) + "," + str(P[2]) + ")")
print("\n")




P_x3_given_x1_x2 = norm.pdf(x = x[3] ,loc = Mu[2], scale = math.sqrt(P_t_tmin1[3] + sigma))
ln_P_x3_given_x1_x2 = math.log(P_x3_given_x1_x2)
print("Answer Q3 ================================")
print("P(z3|x1,x2)=")
print("N(z3;" + str(Mu[2]) + "," + str(P_t_tmin1[3]) + ")")
print("lnP(x3|x1,x2)=")
print(ln_P_x3_given_x1_x2)
print("lnP(x1,x2,x3)=")
ln_P_x1_x2_x3 = ln_P_x3_given_x1_x2 + ln_P_x1_x2
print(ln_P_x1_x2_x3)
print("P(z3|x1,x2,x3)=")
print("N(z3;" + str(Mu[3]) + "," + str(P[3]) + ")")
print("\n")


P_x4_given_x1_x2_x3 = norm.pdf(x = x[4] ,loc = Mu[3], scale = math.sqrt(P_t_tmin1[4] + sigma))
ln_P_x4_given_x1_x2_x3 = math.log(P_x4_given_x1_x2_x3)
print("Answer Q4 ================================")
print("P(z4|x1,x2,x3)=")
print("N(z4;" + str(Mu[3]) + "," + str(P_t_tmin1[4]) + ")")
print("lnP(x4|x1,x2,3)=")
print(ln_P_x4_given_x1_x2_x3)
print("lnP(x1,x2,x3)=")
ln_P_x1_x2_x3_x4 = ln_P_x4_given_x1_x2_x3 + ln_P_x1_x2_x3
print(ln_P_x1_x2_x3_x4)
print("P(z4|x1,x2,x3,x4)=")
print("N(z4;" + str(Mu[4]) + "," + str(P[4]) + ")")
print("\n")




#
# グラフ描画
# if Mac, Duplicated!
#
import matplotlib.pyplot as plt
print("Answer Q5 ================================")
# observed
plt.subplot(3, 1, 1)
plt.plot([x[1],x[2],x[3],x[4]])
# plt.show()

# before update
plt.subplot(3, 1, 2)

print("before updated")
# predicted z1
err1 = 0
print(Mu[0] - err1)
print(Mu[0])
print(Mu[0] + err1)
print("\n")
# predicted z2
err2 = math.sqrt(P_t_tmin1[2]) / 2.0
print(Mu[1] - err2)
print(Mu[1])
print(Mu[1] + err2)
print("\n")
# predicted z3
err3 = math.sqrt(P_t_tmin1[3]) / 2.0
print(Mu[2] - err3)
print(Mu[2])
print(Mu[2] + err3)
print("\n")
# predicted z4
err4 = math.sqrt(P_t_tmin1[4]) / 2.0
print(Mu[3] - err4)
print(Mu[3])
print(Mu[3] + err4)
print("\n")

plt.plot([Mu[0],Mu[1],Mu[2],Mu[3]])
plt.plot([Mu[0]-err1,Mu[1]-2,Mu[2]-err3,Mu[3]-err4])
plt.plot([Mu[0]+err1,Mu[1]+2,Mu[2]+err3,Mu[3]+err4])
# plt.show()

# after update
plt.subplot(3, 1, 3)

print("after updated")
# updated z1
err1 = math.sqrt(P[1]) / 2.0
print(Mu[1] - err1)
print(Mu[1])
print(Mu[1] + err1)
print("\n")
# updated z2
err2 = math.sqrt(P[2]) / 2.0
print(Mu[2] - err2)
print(Mu[2])
print(Mu[2] + err2)
print("\n")
# updated z3
err3 = math.sqrt(P[3]) / 2.0
print(Mu[3] - err3)
print(Mu[3])
print(Mu[3] + err3)
print("\n")
# updated z4
err4 = math.sqrt(P[4]) / 2.0
print(Mu[4] - err4)
print(Mu[4])
print(Mu[4] + err4)
print("\n")


plt.plot([Mu[1],Mu[2],Mu[3],Mu[4]])
plt.plot([Mu[1]-err1,Mu[2]-2,Mu[3]-err3,Mu[4]-err4])
plt.plot([Mu[1]+err1,Mu[2]+2,Mu[3]+err3,Mu[4]+err4])


plt.show()