Implicit particle filters for data assimilation (Q616006)

The authors deal with a model of the state of a system which consists of an Itô stochastic differential equation (SDE): \(dx = f(x,t)dt + g(x,t)dw\), where \(x = (x_1, x_2,\dots,x_m)\) is an \(m\)-dimensional vector, \(f(x,t)\) is an \(m\)-dimensional vector function, \(g(x,t)\) is an \(m\times m\) matrix, and \(w\) is a Brownian motion which encapsulates all the uncertainty in the model. The initial state \(x(0)\) is given and may be random as well. The SDE is supplemented by measurements \(b^n\) at times \(t^n,n = 0,1,\dots\), which are related to the state \(x(t)\) by the relation \(b^n = h(x^n) + GW^n\), where \(h\) is a \(k\)-dimensional, generally nonlinear, vector function with \(k <m\), \(G\) is a matrix, \(x^n = x(t^n)\), and \(W^n\) is a vector whose components are i.i.d. Gaussian variables of mean 0 and variance 1. The filtering problem is to assimilate the data, that is, estimate \(x\) on the basis of both equation and the observations. The authors first discretize the SDE by a difference scheme, so that the equation becomes a discrete recurrence, and assume that the time step in the dynamics equals the fixed time \(\delta\) between observations. For illustration they use the Euler scheme \(x^{n+1} = x^n + \delta f(x^n,n\delta)+ V\), where \(V\) is a Gaussian of mean zero and variance \(g^2(x^n, n\delta)\delta\). Higher-order schemes are also considered. The conditional probability densities \(P(x^n)\) at times \(t^n\), determined by the discretized SDE given the observations, satisfy the recurrence relation \(P(x^{n+1}) = P(x^n)P(x^{n+1}|x^n)P(b^{n+1}|x^{n+1})/Z\), where \(P(x^n) = P(x^n|b^1,b^2,\dots,b^n)\) is the probability density at time \(n\delta\) given the data up to time \(n\delta\), \(P(x^n)P(x^{n+1}|x^n)\) is the probability density of \(x^{n+1}\) given \(x^n\) as it is determined by the dynamics, \(P(b^{n+1}|x^{n+1})\) is the probability of the next observation given the new position, as per the observation equation, and \(Z\) is a normalization constant. The authors estimate \(P_{n+1}\) with the help of \(M\) particles, with positions \(X_i^n\) at time \(t^n\) and \(X_i^{n+1}\) at time \(t^{n+1} (i =1\dots,M)\), which define empirical densities \(\hat{P}_{n}\), \(\hat{P}_{n+1}\) that approximate \(P_{n}\), \(P_{n+1}\) by requiring that, when a particle moves from \(X_i^{n}\) to \(X_i^{n+1}\), the probability of \(X_i^{n+1}\) given \(b^k\) for \(k\leq{n+1}\) be given by \(P(X_i^{n+1}) = P(X_i^n)P(X_i^{n+1}|X_i^n)P(b^{n+1}|X_i^{n+1})/Z_0\), where the hats are omitted, \(P(X_i^n)\), the probability of \(X_i^n\) given \(b^k\) for \(k\leq{n}\), is assumed given, \(P(X_i^{n+1}|X_i^n)\) is determined by the discretized SDE, \(P(b^{n+1}|X_i^{n+1})\) is determined by the observation equation, and \(Z_0\) is an unknown normalization constant. This equation defines the probability density function (pdf) needed to sample for each particle. The idea in implicit sampling is to define probabilities first, and then look for particles that assume them. This is done by choosing a fixed random variable with a given pdf, say a Gaussian \(\exp(-\xi^{\top}\xi/2)/(2\pi)^{m/2}\), which one knows how to sample, and then making \(X_i^{n+1}\) a function of \(\xi\) different for each particle and each step, each of these functions designed so that the map \(\xi\to X_i^{n+1}\) connects highly probable values of \(f\) to highly probable values of \(X_i^{n+1})\). Then, solve the equation \(F_i(X)-\phi_i=\xi^{\top}\xi\), where \(\xi\) is a sample of the fixed variable and \(\phi_i\) is factor needed to make the equation solvable. Once the function \(X_i=X_i(\xi)\) is determined, each value of \(X_i^{n+1}\) appears with probability \(\exp(-\xi^{\top}\xi/2)J^{-1}/(2\pi)^{m/2}\), where \(J\) is the Jacobian of the map \(X_i=X_i(\xi)\), while \(P(X_i^{n+1}|X_i^n)P(b^{n+1}|X_i^{n+1})\) evaluated at \(X_i^{n+1}\) equals \(\exp(-\xi^{\top}\xi/2)\exp{-\phi}\). The sampling weight for the particle is \(\exp{-\phi}J/(2\pi)^{m/2}\). This construction makes it possible to guide the particles to the high-probability area one by one under a wide variety of cases. The equation that links \(\xi\) to \(X\) is underdetermined and its solution can be adopted for each particular problem. The authors present a detailed description of these filters, with illustrative examples, together with new, more general, methods for solving algebraic equations and with a new algorithm for parameter identification.