I need some hints for solving Ecercise 4.4 from Bayesian Filtering & Smoothing by Simo Särkkä:
Select a finite interval in the state space, say, $x \in [-10, 10]$ and discretize it evenly to N subintervals (e.g. N = 1000). Using a suitable numerical approximation to the integrals in the Bayesian filtering equations, implement a finite grid approximation to the Bayesian filter for the Gaussian random walk in Example 4.1. Verify that the result is practically the same as that of the Kalman filter above.
The solution I envisage is:
start from an initial distribution for the hidden state over the intervals;
propagate that discrete distribution and update it such that at each step it incorporates new information.
In a way, just as a particle filter without resampling, only with weight updating, right?
My understanding is that I do not need to use the Kalman filter equations, all I am concerned with is the Bayesian filtering equations of section 4.2.
Here is my code:
import numpy as np from scipy import stats import plotly.graph_objects as go N = 100 xs = np.arange(N) np.random.seed(3) # Simulate Random Walk X = np.zeros(N) for i in range(1, N): X[i] = X[i-1] + np.random.normal() Y = X + np.random.normal(size=N) # Discretization x = np.linspace(-20,20,1001) _x_ = 0.5*(x[1:]+x[:-1]) # use mid-interval values # Initial weight and assumed variance w = np.full(len(_x_), fill_value=1/len(_x_)) P = 2 MM = np.zeros(len(Y)) for i in range(len(Y)): CDF = stats.norm.cdf(x, loc=Y[i], scale=np.sqrt(P+1)) w *= np.diff(CDF) # update distribution w /= np.sum(w) # normalize # w += 1e-16 # avoid degeneracy m = w @ _x_ P = w @ ((_x_-m)**2) MM[i] = m fig = go.Figure() fig.add_scatter(x=xs, y=X, name='Hidden State') fig.add_scatter(x=xs, y=Y, mode='markers', name='Measurements') fig.add_scatter(x=xs, y=MM, name='Inferred') fig.layout.update(height=600) fig.show()
The resulting plot is:
This is clearly not "practically the same as that of the Kalman filter" (provided solution on page 59).
Slightly better, but still, not quite the same.