particle filters with variable number of sampling/factors in models

[iffsid]

A model that involves incremental sampling and factoring of elements, where the number of elements and factors itself can vary across particles, results in the construction of a marginal distribution that has undefined as an element of its support.

The example below illustrates this problem. When the variable numl in inner is set to a constant, the the code runs as expected. However, replacing the number by a call to a stochastic function that varies the number of elements produced, causes the marginal distribution have undefined in its support.

var ds = [0,1,2,3,4,5,6,7,8,9]

var constructL = function(n,m,_l,pf) {
  var _l = _l == undefined ? [] : _l
  var pf = pf == undefined ?  0 : pf
  if (n == 0) {
    factor(-pf)
    return _l
  } else {
    var e = ds[randomInteger(ds.length)]
    var nl = _l.concat([e])
    var nf = m(nl) ? 0 : -30
    factor(nf - pf)
    return constructL(n-1,m,nl,nf)
  }
}

var inner = function(m){
  ParticleFilter(function(){
    var numl = 4 // binomial(0.7,4)+1
    var nl = constructL(numl,m)
    factor(m(nl) ? 0 : -Infinity)
    return nl
  }, 10)
}

var evenM = function(l) {
  // checks if every element in the list is even
  var f = function(a,b){ return a & b }
  var g = function(v) {return v%2 == 0}
  return reduce(function(a,b) { return f(g(a),b) }, g(l[l.length-1]), l.slice(0,-1))
}
var d = inner(evenM)
console.log(d.support([]))
print(d)

[stuhlmueller]

I'm reopening this for now until the issues mentioned in this comment are resolved.

[null-a]

Was this closed by 8688ffa9? (We have #118 for the PF with rejuvenation case.)

[iffsid]

Not entirely. The commits referenced deal with one particular case of variable-number-of-factors; one where the particles are all in sync across observations except at the end, where some of the particles can have extra observations. The more general case of what to do when, say, one particle takes off on some sequence of observations different from other particles, and then syncs back at some later observation, is still not dealt with.

[stuhlmueller]

Aren't all cases of variable numbers of factors such that you have extra factors at the end for some particles? Yes, you might synchronize factors that have very different meanings, which could lead to bad statistical efficiency, but is there a problem with correctness?

[null-a]

I'm trying to understand how we handle this so I can make sure I get it right in #64. If anyone can answer either of the following I'd find it helpful:

• Ignoring rejuvenation, can we say for which models we think the current implementation is correct? (For both the distribution and the normalizing constant estimate.) There seems to be confusion about this in the previous two comments. I'm starting to think it might be correct in general but I only have a rough sketch of an argument which needs more thought.
• I'm not so confident that rejuvenating particles which are at different factor statements is correct. In the case where all particles are at the same factor statement, and we rejuvenate with a kernel which is the program truncated at that factor statement, it makes sense that we end up with particles from the distribution we're interested in. This seems to break down though when we have particles at different factor statements. Is this why we don't test varFactors1 with rejuvenation?

Thanks!

[stuhlmueller]

No rejuvenation

My guess is that the current implementation (in particlefilter.js) results in the correct distribution for any model.

Consider running the particle filter without intermediate resampling steps, only doing a single resampling step at the end based on the total weight for each particle. In this case, it's clear that the distribution is (approximately) correct -- our final step simply changes how we represent the distribution, from being represented in terms of (particle counts, factor weights) to being fully represented in terms of particle counts.

What happens if you introduce an intermediate resampling step? We can think of the original representation as (particle counts, [weights1 + weights2]); the first resampling step changes the representation of the distribution such that the information in weights1 is now represented using particle counts. The second resampling step includes the information in weights2 in the counts representation as well. This is a general argument that intermediate resampling steps are okay.

For this argument, it doesn't matter how many factors (if any) resulted in the weights, whether the number of factors differs between particles, or whether we already applied resampling steps previously. We can always apply resampling to transform the part of the distribution that's represented using weights to being represented using particle counts.

I don't have very high confidence in this argument (both in its correctness and in the fact that it applies to our implementation), but I think it's probably right. I'd guess that a similar argument applies to the normalization constant estimate, but I haven't thought much about it and so have even less confidence there.

Rejuvenation

A first reaction might be: "A particle is approximately sampled from some distribution, and all rejuvenation does is increase the quality of this sample, i.e., make it less approximate. How could this not be correct?"

For models with structural change, this is subtle, though. In our current implementation, the distribution that the particle is sampled from ("program truncated at the factor statement") isn't well-defined for some models.

Consider a model like this:

if (flip(.5)){
  factor(-1);
} else {
  factor(-2);
}
// ... more program here ...

Suppose we run our particle filter and get a particle that arrives at factor(-1). When we run rejuvenation on this particle, we may change the value of flip(.5). In this case, we will never make it to factor(-1), our test for when to stop rejuvenation never succeeds, and we run the program till the end. If we use the trace that results from such rejuvenation in our particle filter, we essentially skip over some factors (we don't accumulate their weights), so I don't expect the distribution to be correct.

One way around this problem would be to only allow proposals to non-structural choices in rejuvenation steps (assuming we can detect such choices, e.g. using static analysis).

[null-a]

so I don't expect the distribution to be correct. One way around this problem would be to only allow proposals to non-structural choices in rejuvenation steps

To alert users to the possible incorrectness in the interim, should we consider throwing an error (or warning) if we try to rejuvenate particles which are out of sync? (i.e. at different factor statements.)

[stuhlmueller]

Yes, I think it's a good idea to throw an error in this case.

[iffsid]

Does this mean it will now throw an error in the extra-factors-at-the-end situation that we currently handle? If so, could we have the option to override the error throwing? I'm thinking the example @andreas mentions as being suspect in the rejuvenation case could potentially be handled by using setProposalBoundary.

// ....
var ifTest = flip(0.5)
setProposalBoundary();
if (ifTest){
  // some code
  factor(-1);
} else {
  // some other code
  factor(-2);
}
// ....

[null-a]

Does this mean it will now throw an error in the extra-factors-at-the-end situation

By extra-factors-at-the-end do you mean something like extra factor statements at the end of the program without intervening sample statements? For example:

var x = flip(0.5);
factor(-1);
if (x) {
  factor(-2);
}
// End of program

that we currently handle?

Do you mean we're confident that rejuvenation in this case doesn't run into the problem Andreas mentioned?

[iffsid]

Well there could be intervening sample statements

var x = flip(0.5);
setProposalBoundary()
var a1 = flip(0.2);
factor(a1 ? 0 : -Infinity);
if (x) {
  var a2 = flip(0.4);
  factor(a2 ? 0 : -Infinity);
}

The point being that if you could move the structural random choices outside the rejuvenation boundary, then the fact that you could have particles with one or two factors wouldn't break things.

I wanted the option to override throwing the error because there are cases (as above, and when constructing sequences involving heuristic factors) where rejuvenating particles that appear out of sync is not wrong.

[stuhlmueller]

I'm not opposed to having an option (such as allowDangerousOutOfSyncRejuvenation) that disables throwing the error. This would be temporary and would be removed once we implement a coherent strategy for rejuvenation in models with structural change.

[null-a]

One way around this problem would be to only allow proposals to non-structural choices in rejuvenation steps

I'm not sure that's enough, assuming the rest of the algorithm is unchanged. Consider this model:

var model = function() {
  var x = flip(0.5);
  var y = flip(0.5);
  setProposalBoundary()
  var z = flip(0.5);
  if (x) {
    factor(x ? 0 : -1);
  }
  if (y) {
    factor(y ? 0 : -1);
  }
  return [x, y]
};

Here, I think we still have particles failing to exit at the correct limitAddress during rejuvenation, even though we don't propose to x or y. This is because not all particles see every factor statement but we currently use the same limitAddress for each particle during rejuvenation.

If that's correct, it suggests setting limitAddress per-particle to the address of the last factor statement reached.

(I'm thinking about this because in #64 I assert that particles reach limitAddress during rejuvenation, and I'm trying to understand whether I need to relax that.)

[stuhlmueller]

If that's correct, it suggests setting limitAddress per-particle to the address of the last factor statement reached.

Yes, if we want to allow use of SMC for models with structural change (even if the actual rejuvenation steps only propose to non-structural choices), we would have to keep track of limitAddress per particle.

(Just to be clear: Since we don't have a way of limiting proposals to non-structural choices for now, we should continue throwing an error by default if we try to rejuvenate particles that are out of sync. Working out rejuvenation for models with structural change is a separate task from getting plain particle filtering working - now issue #176.)

[stuhlmueller]

If that's correct, it suggests setting limitAddress per-particle to the address of the last factor statement reached.

As Sid just pointed out to me: the current implementation does this.

[null-a]

As Sid just pointed out to me: the current implementation does this.

Oops, sorry! I did check before suggesting it, I must have misunderstood what's happening.

I looked here and thought that it looked like rejuvenation was using limitAddress = a where a is the address of the factor statement reached by the final particle. What am I missing?

[stuhlmueller]

No, sorry, you're right. I should have checked -- turns out it's just Sid's local branch that saves the current factor address for each particle and then use that as limit address.

[null-a]

we should continue throwing an error by default if we try to rejuvenate particles that are out of sync

I've added an allowOutOfSyncRejuv option to the SMC implementation in #169 that does this. It also now tracks a limitAddress per-particle. And just to clarify, it's still the case that it checks that particles reach their limitAddress during rejuvenation, and errors otherwise.

I don't have any more related changes planned for #169. @iffsid Are you happy this leaves things in a usable state for you?

[iffsid]

Yeah, looks like that will help. Thanks! I'll take a closer look at things soon; haven't followed the changes over the last couple of weeks.

[null-a]

I've added an allowOutOfSyncRejuv option to the SMC implementation

To clarify, as of #245 this option no longer exists and "out of sync" rejuvenation is allowed by default.

Are we ok to close this? (If not, we should maybe open new tickets for specific issues as this thread is pretty convoluted at this point.)