Example 5.8#
Source localization example.
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
rng = np.random.default_rng(32)
x_true = np.array([0, 1])
R = 0.1
s = np.array([[2, 1.2], [-1, -1],[-0.5, 3]])
def p(x, mu, K): # multivariate gaussian density using scipy stats
return stats.multivariate_normal.pdf(x, mean=mu, cov=K)
def log_p(x, mu, K): # log of multivariate gaussian density using scipy stats
return stats.multivariate_normal.logpdf(x, mean=mu, cov=K)
def p_1d(x, mu, sig): # 1d gaussian
return 1/np.sqrt(2*np.pi * sig**2) * np.exp(-(x-mu)**2/(2*sig**2))
def log_p_1d(x, mu, sig):
return -0.5*np.log(2*np.pi*sig**2) - (x-mu)**2/(2*sig**2)
y = np.zeros((3, 1))
y[0] = np.random.normal(np.linalg.norm(s[0] - x_true), R, 1)
y[1] = np.random.normal(np.linalg.norm(s[1] - x_true), R, 1)
y[2] = np.random.normal(np.linalg.norm(s[2] - x_true), R, 1)
N = 5000
m = np.array([0, 0])
P = 20 * np.array([[1, 0], [0, 1]])
sigma = 0.04
P_prop = sigma * np.array([[1, 0], [0, 1]])
x = np.zeros((2, N))
x[:, 0] = np.array([-3, -3])
burnin = 500
acc = 0
fig = plt.figure(figsize=(15, 7))
for n in range(1, N):
x_s = rng.multivariate_normal(x[:, n-1], P_prop, 1)
u = rng.uniform(0, 1, 1)
# we will terms separately
logr1 = log_p(x_s, m, P) - log_p(x[:, n-1], m, P) # terms relating to p(x)
logr2 = 0
for i in range(3):
# compute the lik term relating to the datapoint i
logr2_i = log_p_1d(y[i], np.linalg.norm(s[i] - x_s), R) - log_p_1d(y[i], np.linalg.norm(s[i] - x[:, n-1]), R)
logr2 = logr2 + logr2_i
logr = logr1 + logr2
if np.log(u) < logr:
x[:, n] = x_s
acc += 1
else:
x[:, n] = x[:, n-1]
plt.clf()
plt.subplot(1,2,1)
plt.scatter(x[0, burnin:], x[1, burnin:], color='blue', alpha=0.5, label='Samples')
plt.scatter(x_true[0], x_true[1], color='k', marker='o', label='True')
plt.scatter(s[0][0], s[0][1], color=[0.8, 0, 0], marker='o', label='Sensor')
plt.scatter(s[1][0], s[1][1], color=[0.8, 0, 0], marker='o')
plt.scatter(s[2][0], s[2][1], color=[0.8, 0, 0], marker='o')
plt.plot([x_true[0], s[0][0]], [x_true[1], s[0][1]], color='k', linestyle='--')
plt.plot([x_true[0], s[1][0]], [x_true[1], s[1][1]], color='k', linestyle='--')
plt.plot([x_true[0], s[2][0]], [x_true[1], s[2][1]], color='k', linestyle='--')
plt.legend()
plt.title('The object, the sensors, the posterior')
plt.subplot(1,2,2)
plt.plot(x[0, burnin:], color=[0.8, 0, 0], alpha=0.5)
plt.plot(x[1, burnin:], '.--', color=[0.8, 0, 0], markersize=0.1, alpha = 0.5)
plt.title('Path traces of the 2D Markov chain')
plt.show()
