Commit 0f6f0a8b authored by Guillaume Chérel's avatar Guillaume Chérel
Browse files

ABC: Avoid generating points with null prior.

parent 3ddd7b0e
......@@ -107,7 +107,7 @@ object APMC {
def exposedStep(p: Params)(implicit rng: util.Random): ExposedEval[State, Matrix, (State, Matrix, Matrix), Matrix, State] =
ExposedEval(
pre = { s =>
val (sigmaSquared, newThetas) = stepPreEval(p.n, p.nAlpha, s)
val (sigmaSquared, newThetas) = stepPreEval(p.n, p.nAlpha, p.priorDensity, s)
((s, sigmaSquared, newThetas), newThetas)
},
post = { (sstep, newXs) =>
......@@ -115,7 +115,11 @@ object APMC {
stepPostEval(p.n, p.nAlpha, p.priorDensity, p.observed, s, sigmaSquared, newThetas, newXs)
})
def stepPreEval(n: Int, nAlpha: Int, s: State)(implicit rng: util.Random): (Matrix, Matrix) = {
def stepPreEval(
n: Int,
nAlpha: Int,
priorDensity: Array[Double] => Double,
s: State)(implicit rng: util.Random): (Matrix, Matrix) = {
val thetasM = MatrixUtils.createRealMatrix(s.thetas)
val dim = thetasM.getColumnDimension()
......@@ -124,8 +128,10 @@ object APMC {
val weightedDistributionTheta = new EnumeratedIntegerDistribution(apacheRandom(rng), Array.range(0, nAlpha), s.weights)
val newThetas =
Array.fill(n - nAlpha) {
val resampledTheta = thetasM.getRow(weightedDistributionTheta.sample)
new MultivariateNormalDistribution(apacheRandom(rng), resampledTheta, sigmaSquared.getData).sample
Iterator.continually {
val resampledTheta = thetasM.getRow(weightedDistributionTheta.sample)
new MultivariateNormalDistribution(apacheRandom(rng), resampledTheta, sigmaSquared.getData).sample
}.dropWhile { priorDensity(_) == 0 }.next
}
(sigmaSquared.getData, newThetas)
......
......@@ -142,7 +142,12 @@ object MonAPMC {
type StepState = Either[Matrix, (APMC.State, Int, Matrix, Matrix)]
def preStep(n: Int, nAlpha: Int, priorSample: util.Random => Array[Double], state: MonState)(implicit rng: util.Random): (StepState, Matrix) = {
def preStep(
n: Int,
nAlpha: Int,
priorSample: util.Random => Array[Double],
priorDensity: Array[Double] => Double,
state: MonState)(implicit rng: util.Random): (StepState, Matrix) = {
/* If s is empty or the number of particles it contains hasn't reach nAlpha yet, keep generating particles (n - nAlpha by n - nAlpha) from the prior using the function APMC.stepOne and setting the parameter n to (n - nAlpha)*/
val reducedN = n - nAlpha
......@@ -157,13 +162,19 @@ object MonAPMC {
val thetas = APMC.initPreEval(reducedN, priorSample)
(Left(thetas), thetas)
} else {
val (sigmaSquared, thetas) = APMC.stepPreEval(n, nAlpha, s)
val (sigmaSquared, thetas) = APMC.stepPreEval(n, nAlpha, priorDensity, s)
(Right((s, t0, sigmaSquared, thetas)), thetas)
}
}
}
def postStep(n: Int, nAlpha: Int, priorDensity: Array[Double] => Double, observed: Array[Double], stepState: StepState, xs: Matrix)(implicit rng: util.Random) = {
def postStep(
n: Int,
nAlpha: Int,
priorDensity: Array[Double] => Double,
observed: Array[Double],
stepState: StepState,
xs: Matrix)(implicit rng: util.Random) = {
val reducedN = n - nAlpha
stepState match {
case Left(thetas) => State(0, APMC.initPostEval(reducedN, nAlpha, observed, thetas, xs))
......@@ -178,7 +189,7 @@ object MonAPMC {
def exposedStep(p: Params)(implicit rng: util.Random): ExposedEval[MonState, Matrix, Either[Matrix, (APMC.State, Int, Matrix, Matrix)], Matrix, MonState] =
ExposedEval(
pre = preStep(p.apmcP.n, p.apmcP.nAlpha, p.apmcP.priorSample, _),
pre = preStep(p.apmcP.n, p.apmcP.nAlpha, p.apmcP.priorSample, p.apmcP.priorDensity, _),
post = postStep(p.apmcP.n, p.apmcP.nAlpha, p.apmcP.priorDensity, p.apmcP.observed, _, _))
def steps(s: MonState) = s match {
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment