public int compare(Double o1, Double o2) {
return o1 < o2 ? 1 : -1;
}
});
for (Entry<FruehwirthMultiParticle, ? extends Number> particleEntry : target.asMap().entrySet()) {
final FruehwirthMultiParticle particle = particleEntry.getKey();
/*
* Fruewirth-Schnatter's method for sampling...
* TODO do we sample new y's for each particle?
* TODO when using FS aug sampling, we should probably go all the way
* and replicate their beta sampling.
* That would require a change here...
*/
final double U = this.random.nextDouble();
double lambdaSum = 0d;
Vector sampledAugResponse = VectorFactory.getDefault().createVector(this.numCategories);
for (int i = 0; i < this.numCategories; i++) {
// final Vector betaSample = particle.getLinearState().sample(this.random);
final MultivariateGaussian predictivePrior = particle.getLinearState().clone();
KalmanFilter kf = particle.getRegressionFilter(i);
final Matrix G = kf.getModel().getA();
predictivePrior.setMean(G.times(predictivePrior.getMean()));
predictivePrior.setCovariance(
G.times(predictivePrior.getCovariance()).times(G.transpose())
.plus(kf.getModelCovariance()));
// X * beta
final double lambda = Math.exp(data.getObservedData().times(
predictivePrior.getMean()).getElement(0));
lambdaSum += lambda;
final double dSampledAugResponse = -Math.log(
-Math.log(U)/(1d+lambdaSum)
- (data.getObservedValue().getElement(0) > 0d
? 0d : Math.log(this.random.nextDouble())/lambda));
sampledAugResponse.setElement(i, dSampledAugResponse);
}
/*
* Expand particle set over mixture components
*/
for (int j = 0; j < 10; j++) {
double categoriesTotalLogLikelihood = 0d;
final UnivariateGaussian componentDist =
this.evDistribution.getDistributions().get(j);
for (int k = 0; k < this.numCategories; k++) {
/*
* TODO could avoid cloning if we didn't change the measurement covariance,
* but instead used the componentDist explicitly.
*/
final FruehwirthMultiParticle predictiveParticle = particle.clone();
predictiveParticle.setAugResponseSample(sampledAugResponse);
predictiveParticle.setEVcomponent(componentDist);
/*
* Update the observed data for the regression component.
*/
predictiveParticle.getRegressionFilter(k).getModel().setC(data.getObservedData());
final Matrix compVar = MatrixFactory.getDefault().copyArray(
new double[][] {{componentDist.getVariance()}});
predictiveParticle.getRegressionFilter(k).setMeasurementCovariance(compVar);
final double logLikelihood = this.updater.computeLogLikelihood(predictiveParticle, data)
+ Math.log(this.evDistribution.getPriorWeights()[j])
+ (particleEntry.getValue().doubleValue() - prevTotalLogLikelihood);
categoriesTotalLogLikelihood = LogMath.subtract(categoriesTotalLogLikelihood, logLikelihood);
particleTotalLogLikelihood = LogMath.add(particleTotalLogLikelihood, logLikelihood);
particleTree.put(logLikelihood, predictiveParticle);
}
}
}
// final WFCountedDataDistribution<FruehwirthMultiParticle> resampledParticles =
// ExtSamplingUtils.waterFillingResample(logLikelihoods, particleTotalLogLikelihood,
// particleSupport, random, this.numParticles);
final WFCountedDataDistribution<FruehwirthMultiParticle> resampledParticles =
ExtSamplingUtils.waterFillingResample(
Doubles.toArray(particleTree.keySet()),
particleTotalLogLikelihood,
Lists.newArrayList(particleTree.values()),
random, this.numParticles);
/*
* Propagate
*/
target.clear();
for (final Entry<FruehwirthMultiParticle, MutableDouble> particleEntry : resampledParticles.asMap().entrySet()) {
final FruehwirthMultiParticle updatedParticle = sufficientStatUpdate(particleEntry.getKey(), data);
target.set(updatedParticle, particleEntry.getValue().value);
}
Preconditions.checkState(target.getDomainSize() == this.numParticles);
}