Context
For method = "is" / "smc-bpf" / "smc-apf", tp_run()'s raw output per sweep is list(samples = [...], weights = [...], normConst = <scalar>). There's no built-in function that turns this into a tidy "one row per particle, one column per parameter, plus a weight column" data frame the way one would want for:
- computing an unbiased marginal-likelihood (log Z) estimate,
- ESS diagnostics,
- weighted posterior summaries.
tp_parse_smc() exists but (see #19) its weight/norm_constant columns don't correspond to the raw output at all, so it isn't usable for this today.
What we ended up writing
flatten_sample <- function(s) {
d <- s[["__data__"]]
out <- list()
for (nm in names(d)) {
v <- d[[nm]]
if (is.list(v) && length(v) > 1) { # e.g. an array-valued field like `lambda`
lam <- as.numeric(unlist(v))
names(lam) <- paste0(nm, seq_along(lam))
out <- c(out, as.list(lam))
} else {
out[[nm]] <- if (is.list(v)) v[[1]] else v
}
}
out
}
tp_parse_is <- function(out) {
sweeps <- lapply(seq_along(out), function(s) {
samples <- lapply(out[[s]]$samples, flatten_sample)
draws <- do.call(rbind, lapply(samples, as.data.frame))
draws$run <- s
draws$particle <- seq_len(nrow(draws))
draws$logWeight <- as.numeric(unlist(out[[s]]$weights)) # see #20 re: -Inf handling
list(draws = draws, logZ = out[[s]]$normConst)
})
draws <- do.call(rbind, lapply(sweeps, `[[`, "draws"))
logZ_per_run <- vapply(sweeps, `[[`, numeric(1), "logZ")
m <- max(logZ_per_run)
list(draws = draws, logZ = m + log(mean(exp(logZ_per_run - m))))
}
(Combining multiple sweeps' normConst values correctly requires logmeanexp, not a plain average — worth documenting explicitly if this becomes a shipped function, since it's an easy mistake.)
Request
A maintained, correctly-implemented equivalent of the above (fixing #19, or as a new function) would remove a very common pain point for any exact-vs-approximate model validation work, which importance sampling / SMC's normConst is naturally suited for.
Context
For
method = "is"/"smc-bpf"/"smc-apf",tp_run()'s raw output per sweep islist(samples = [...], weights = [...], normConst = <scalar>). There's no built-in function that turns this into a tidy "one row per particle, one column per parameter, plus a weight column" data frame the way one would want for:tp_parse_smc()exists but (see #19) its weight/norm_constant columns don't correspond to the raw output at all, so it isn't usable for this today.What we ended up writing
(Combining multiple sweeps' normConst values correctly requires
logmeanexp, not a plain average — worth documenting explicitly if this becomes a shipped function, since it's an easy mistake.)Request
A maintained, correctly-implemented equivalent of the above (fixing #19, or as a new function) would remove a very common pain point for any exact-vs-approximate model validation work, which importance sampling / SMC's
normConstis naturally suited for.