0
$\begingroup$

So, my endeavor to apply the is just for my own edificationI am currently struggling with an attempt to apply a bootstrap particle filter (Gordon, Salmond, Smith, 2003) to a linear, Gaussian state-space model

$$s_t=A\,s_{t-1}+B\,\nu_t\qquad\text{( transition equation )}$$ $$\qquad z_t=C\,s_t+D\,\varepsilon_t\qquad\text{( measurement equation )}$$

where the state has dimension $n_s$, the state innovations have dimension $n_\nu$, the observables have dimension $n_z$, and the measurement errors have dimension $n_\varepsilon$. Moreover, $\nu_t$ and $\varepsilon_t$ are Gaussian, serially and mutually uncorrelated and have variance-covariance matrix $I$, respectively.

Given the linear and Gaussian nature of my problem, I realize that conventional Kalman filtering methods would be more than sufficient to tackle any problem involving this state-space model. This is just to understand the way the particle filter works.

As far as I understood, the algorithms runs as follows:

I start out with a prior for $s_0$, say $s_{0}\sim\mathcal{N}(0,1)$, and draw a set of random (iid) samples $\{s_{i,0}\}_{i=1}^N$. I sort of pretend that these come from a density $f(s_{0}|z_0)$ and when I advance from $t=1$ to $t=2$ they actually will come (in an approximate sense) from $f(s_1|z_1,z_0).$

Then, I use my transition equation to simulate $s_{i,1}$ for every $i$. That is, I draw a set $\{\nu_{i}\}_{i=1}^N$ from $\mathcal{N}(0,I)$ and evaluate $$s_{i,1}=A\,s_{i,0}+B\,\nu_i$$ which yields a sample $\{s_{i,1}\}_{i=1}^N$.

Given my measurement equation, I know that $z_1|s_{i,1}$ is distributed according to a normal with mean $A\,s_{i,1}$ and variance-covariance matrix $DD^\prime.$ Hence, I can evaluate the pdf $f(z_1|s_{i,1})$ for every single $s_{i,1}$.

Provided that all of the foregoing is correct, this is the point were my problem arises.

In particular, I have $$p(z_1|s_{i,1})=|2\pi DD^\prime|^{-1/2}\exp(\,-1/2\,\cdot\,(z_1-Cs_{i,1})^\prime(DD^\prime)^{-1}(z_1-Cs_{i,1})\,)$$ Given (not unreasonable) assumptions on the dimensionality of the state-space system as well as the "magnitude" of $D$, it is easy to get to a point where the logarithm of conditional density above takes values somewhere in the vicinity of -8000. Hence, when coding the problem in Matlab, the exponential of this is returned as 0. Now, in order to loop over $t$ in this algorithm, I need to obtain an expression for the importance weights$$\dfrac{p(z_1|s_{i,1})}{\sum_{i=1}^Np(z_1|s_{i,1})}$$ such that I cannot simply work with the log-density. How do I proceed, I see that $p(z_1|s_{i,1})$ is not actually zero, but it would seem that I cannot manage to make Matlab see the difference.

I might be at a loss here, but am I missing something? Am I doing something wrong? Is there a trick to that?

Many Thanks.

Edit - 2nd Question

In the following $T$ is the number of data points and $N$ is the number of particles.

My problem is that in every iteration a single particle is assigned an importance weight of 1. While I could get around that by scaling up the variance-covariance matrix of the measurement equation, this does not help with my second problem.

In particular, I am trying to evaluate the likelihood $f(z_{1:T})$ for a given parameterization of $A,B,C,$ and $D$ such that I can run a numerical optimizer to find the parameters that maximize the likelihood. As I simulated the system, I do know the true values.

When I do all of that using the Kalman filter, it is true that for every parameter with the other parameters fixed at their true value) the likelihood attain its peak very close to the true value.

This is not the case, when I use the particle filter.

Assuming even weights (i.e. resetting weights to $1/N$ after resampling), I can approximate $$f(z_{t+1}|z_{1:t})\approx\dfrac{1}{N}\sum_{i=1}^Nf(z_{t+1}|s_{it+1})$$

Since I am only interested in the object $$\log\,f(z_{1:T})=\sum\,\log\,f(z_{t+1}|z_{1:t})$$ I overwrite the filtered states in each iteration of the following piece of code and just update the likelihood

Q=B*B'; % vcv for the innovation of the transition equation
R=D*D'; % vcv for the innovation of the measurement equation

QT=Q';
QT=QT(:); % vectorizing transition equation vcv

Sig0= inv(eye(36)-kron(A,A))*QT(:); 
Sig0= reshape(Sig0,6,6); % uncoditional variance of state

s0=zeros(6,1); % unconditional mean of state

s=mvnrnd(s0,Sig0,N)'; % drawing N samples for initial particle cloud
v=randn(4,N,T); % innovations to propagate particle cloud

L=0; % initializing log-likelihood

for t=1:T  
    % take s from last iteration and propagate via transition equation - returns 6-by-N array 
    s=A*s+B*v(:,:,t);
    % calculate the density p(z_t|s_it) for each s_it in cloud - returns a N-by-1 vector
    for i=1:N
        p(i)=-1/2*[log(det(2*pi*R))+(z(:,t)-C*s(:,i))'*inv(R)*(z(:,t)-C*s(:,i))];
    end
    % obtain max element of p for numerically robust computations 
    pm=max(p);
    % update the log-likelihood such that at the t-th iteration lp is the sum
    % from k=1:t of the means (over i) of p(z_k|s_ik).
    % Since I take an unweighted mean, this is as if I would have made
    % weights even after resampling in the last step. 
    L=L-log(N)+pm+log(sum(exp(p-pm)));
    % calculate the importance weight
    w=exp(p-pm)./sum(exp(p-pm));
    % multinomial resampling at each step according to importance weights
    id=randsample(1:N,N,true,w);
    x_s=x_s(:,id);
end

I would appreciate any pointer that might help me figuring out what is wrong with my code.

Many thanks.

$\endgroup$
  • $\begingroup$ What is the dimensionality of your state space? Particle filters do not scale well to high-dimensional problems, in the first place. Nonetheless, you appear to be doing bootstrap filtering, in which case I don't know why A is in your weight function $p(z_1|s_{i,1})$. $\endgroup$ – Forgottenscience Oct 31 '19 at 18:58
  • $\begingroup$ Thank you very much for your message. For the purpose of my endeavor to gain a bit of facility with the particle filter, I simulated a system with a 6-dimensional state-space. Regarding the $A$ in my weight function, that is just a typo. Sorry about that. I edited my post accordingly. $\endgroup$ – YL-Wint Oct 31 '19 at 19:03
  • $\begingroup$ It happens all the time. Use the log sumexp trick $\endgroup$ – Taylor Oct 31 '19 at 19:31
  • $\begingroup$ stats.stackexchange.com/questions/304758/softmax-overflow/… $\endgroup$ – Taylor Oct 31 '19 at 19:40
  • 1
    $\begingroup$ Thank you very much. That is the trick, I was looking for. $\endgroup$ – YL-Wint Oct 31 '19 at 20:11
1
$\begingroup$

For simplicity and to follow the classical state space filtering literature, I will denote by $f(s_t|s_{t-1})$ the transition density of the state $s_t$ given $s_{t-1}$ and $g(z_t|s_t)$ the measurement density of the observation $z_t$ given $s_t$. Independent of your transition and measurement equations, to run the vanilla bootstrap particle filter, you do the following:

We have $N$ particles, $T$ time steps.

At the "initial time" $t=0$:

  1. Draw initial samples $s_0^n \sim p(s_0^n)$ for $n = 1,2,\ldots, N$, we denote the collection of particles by $(s_0^n)$.

for $t = 1,2, \ldots T$:

  1. Propagate your sample cloud $(s_{t-1}^n)$ to time $t$ by the transition equation.
  2. For each $n=1,2,\ldots,N$ in your cloud calculate $w_t^n = g(z_t|s_t^n)$, these are your particle weights.
  3. Normalize the weights $W^n_t = \frac{w^n_t}{\sum_n w^n_t}$
  4. Resample your particles with probability of drawing the $i$'th particle proportional to its weight. In practice this is $N$ samples, denote the $n$'th of them by $\phi_t^n$, from a categorical distribution with the probability of choosing index $i$ given by $W^i_t$, that is, $P(\phi_t^n = i) = W^n_t$.
  5. Replace your point cloud $(s_t^n)$ with $(s_t^{\phi_t^n})$ and set each weight equal to $\frac 1 N$.

Given the dimensionality of your example and that it is purely Gaussian, the particle filter will most likely be doing a very good job if implemented correctly. My guess is that you have a coding error somewhere when calculating the weights. If you ever try to run a particle filter in high dimensions and literally cannot deal with the exponential probabilities of the weights, you can use the exp-normalize trick.

$\endgroup$
  • $\begingroup$ Thank you so very much for your answer. Makes perfect sense. I tried to implement my algorithm according to the steps you delineated above and it would seem that in every iteration, I have a single importance weight equal to one. So I suppose your guess that I might have a coding error somewhere in there is correct. I'll dive into that. $\endgroup$ – YL-Wint Oct 31 '19 at 20:10
  • $\begingroup$ If I may ask another (related) question. When I initialize my procedure, I do so using a draw from a normal with the unconditional mean and vcv of $s_0$. Then, I compute my particle weight where, upon exp-normalization, it turns out that a single particle is assigned a weight of 1. Hence, my resampled particle cloud consists of $N$ copies of that particle. When I propagate this cloud and compute the corresponding weights, once again, one particle has weight 1. Does that sound like a coding error? Or could that be due to the parameterization of my problem (which is calibrated to fit real data)? $\endgroup$ – YL-Wint Nov 1 '19 at 2:54
  • $\begingroup$ While PF's do collapse in high dimensions, you should not be able to have that poor coverage in 6 dimensions. Even if you only had $N=2$ I would imagine you'd be hard-pressed to have a single particle take all the relative weight. It must be an error somewhere. $\endgroup$ – Forgottenscience Nov 1 '19 at 10:25
  • $\begingroup$ Since I am at my wits end, I edited my question and posted the segment of my code containing the error. I tried to annotate my code and explain my overall goal as concisely as I could. I'd very much appreciate, I you could have a look. Thank you. $\endgroup$ – YL-Wint Nov 1 '19 at 23:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.