turinglang / advancedps.jl Goto Github PK

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Following the discussion in TuringLang/Libtask.jl#142 we need to update the split logic slightly since the reference are update in Libtask.tape_copy when we fork tasks.

forkr is likely broken without explicit dependence on the replaying mechanism in Turing.jl. This is required to correctly implement the reference particle for particle Gibbs. A potential fix is the switch to a simpler replying mechanism by fixing the random seed for the reference particle.

Does this package support multidimensional state-spaces?

Reading the code it doesn't seem so, and trying to run the following block with a 2D model it doesn't seem so, but I'm new to Julia programming and unsure whether I have bugs.

model = Log2DRatePoissonSSM(θ₀)
pgas = AdvancedPS.PGAS(Nₚ)
chains = sample(rng, model, pgas, Nₛ; progress=false);

As originally highlighted in the discussion reported below, I am trying to use Turing.jl for online learning. I have noticed that the SMC's implementation is based on AdvancedPS.jl and - given that the relevant part of the documentation is missing - decided to open an issue to have the following clarifications.

1. Suppose that I have estimated a Turing model via SMC by calling the sample(...) function and on the basis of a univariate time series with T datapoints. Assume now that I have additional data for the following h points in time and that I would like to update my estimation. How should I do it starting from the sample(...) output?

2. How do I specify a custom importance distribution from Turing?

3. Any detail specific to the example mentioned below.

   Hi there, I am trying to use Turing.jl for online learning (using SMC). For now, I am trying to replicate the example in Section 5.2 (Mixtures) of Nicolas Chopin (2002 - link here) to have a better understanding of Turing's syntax and scope. I cannot find a dedicated example for online learning in the documentation. Is there any external tutorial / guide I should read? If not, how would you suggest to implement an example similar to the one mentioned above?

Originally posted by @fipelle in TuringLang/Turing.jl#1895

Working through this example, and it suggests using PGAS to help with sample impoverishment. I've tried only replacing pgas = AdvancedPS.PG(n_particles) with pgas = AdvancedPS.PGAS(n_particles) as is done here , while keeping all else the same but am met with this error

ERROR: MethodError: no method matching step(::MersenneTwister, ::NonLinearTimeSeries, ::AdvancedPS.PGAS{AdvancedPS.ResampleWithESSThreshold{typeof(resample_systematic), Float64}})
Closest candidates are:
  step(::AbstractRNG, ::AdvancedPS.AbstractStateSpaceModel, ::AdvancedPS.PGAS) at ~/.julia/packages/AdvancedPS/Vox9w/src/smc.jl:130
  step(::AbstractRNG, ::AdvancedPS.AbstractStateSpaceModel, ::AdvancedPS.PGAS, ::Union{Nothing, AdvancedPS.PGState}; kwargs...) at ~/.julia/packages/AdvancedPS/Vox9w/src/smc.jl:130
  step(::AbstractRNG, ::AbstractMCMC.AbstractModel, ::AdvancedPS.PG) at ~/.julia/packages/AdvancedPS/Vox9w/src/smc.jl:88
  ...
Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/AbstractMCMC/fWWW0/src/sample.jl:125 [inlined]
  [2] macro expansion
    @ ~/.julia/packages/ProgressLogging/6KXlp/src/ProgressLogging.jl:328 [inlined]
  [3] (::AbstractMCMC.var"#21#22"{Bool, String, Nothing, Int64, Int64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, MersenneTwister, NonLinearTimeSeries, AdvancedPS.PGAS{AdvancedPS.ResampleWithESSThreshold{typeof(resample_systematic), Float64}}, Int64, Int64})()
    @ AbstractMCMC ~/.julia/packages/AbstractMCMC/fWWW0/src/logging.jl:12
  [4] with_logstate(f::Function, logstate::Any)
    @ Base.CoreLogging ./logging.jl:511
  [5] with_logger(f::Function, logger::LoggingExtras.TeeLogger{Tuple{LoggingExtras.EarlyFilteredLogger{TerminalLoggers.TerminalLogger, AbstractMCMC.var"#1#3"{Module}}, LoggingExtras.EarlyFilteredLogger{Logging.ConsoleLogger, AbstractMCMC.var"#2#4"{Module}}}})
    @ Base.CoreLogging ./logging.jl:623
  [6] with_progresslogger(f::Function, _module::Module, logger::Logging.ConsoleLogger)
    @ AbstractMCMC ~/.julia/packages/AbstractMCMC/fWWW0/src/logging.jl:36
  [7] macro expansion
    @ ~/.julia/packages/AbstractMCMC/fWWW0/src/logging.jl:11 [inlined]
  [8] mcmcsample(rng::MersenneTwister, model::NonLinearTimeSeries, sampler::AdvancedPS.PGAS{AdvancedPS.ResampleWithESSThreshold{typeof(resample_systematic), Float64}}, N::Int64; progress::Bool, progressname::String, callback::Nothing, discard_initial::Int64, thinning::Int64, chain_type::Type, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ AbstractMCMC ~/.julia/packages/AbstractMCMC/fWWW0/src/sample.jl:116
  [9] sample(rng::MersenneTwister, model::NonLinearTimeSeries, sampler::AdvancedPS.PGAS{AdvancedPS.ResampleWithESSThreshold{typeof(resample_systematic), Float64}}, N_or_isdone::Int64; kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:progress,), Tuple{Bool}}})
    @ AbstractMCMC ~/.julia/packages/AbstractMCMC/fWWW0/src/sample.jl:51
 [10] top-level scope
    @ REPL[35]:1

could someone help with how to resolve this?

See TuringLang/Turing.jl#1858 and #69

Seems like the latest release of Libtask broke something here:
https://github.com/TuringLang/AdvancedPS.jl/runs/6924872237?check_suite_focus=true

See here. This will be fixed in new versions of Libtask - hopefully soon. We should warn the user when making a new Turing release - examples depend on this feature will return incorrect results but won't throw errors.

Also, unbounded loops are not supported but hopefully will be fixed in Libtask soon.

I'm wondering if a model like the one I present below is possible? The basic problem is one where the state isn't a single distribution, but a collection of distributions, which all evolve in a Markovian manner. I don't know exactly how this works internally, so the code is based on the assumption that the state is propagated forward from initialization to transition to observation.

n_trials, n_cols = size(X)

Parameters = @NamedTuple begin
    X::Matrix
    lam_lapse_init::Float64
    sigma_set_init::Array{Float64}
    mu_init::Array{Float64}
    n_trials::Int64
    n_cols::Int64
end

mutable struct PF <: AdvancedPS.AbstractStateSpaceModel
    W::Matrix
    lam_lapse::Array
    sigma_set::Matrix
    theta::Parameters
    PF(theta::Parameters) = new(
        zeros(Float64, theta.n_trials, theta.n_cols),
        zeros(Float64, theta.n_trials),
        zeros(Float64, theta.n_trials, theta.n_cols),
        theta
    )
end


function init_step(m::PF)
    return (
        truncated(Normal(m.theta.lam_lapse_init, 0.1), lower=-10),
        truncated(Normal(m.theta.sigma_set_init[1], 0.1), lower=0.),
        truncated(Normal(m.theta.sigma_set_init[2], 0.1), lower=0.),
        truncated(Normal(m.theta.sigma_set_init[3], 0.1), lower=0.),
        MvNormal(m.theta.mu_init, 1.)
    )
end

AdvancedPS.initialization(m::PF) = init_step(m)

function transition_step(m::PF, state)
    return (
        truncated(Normal(state[1], 0.1), lower=-10), # lam_lapse
        truncated(Normal(state[2], 0.1), lower=0.), # sigma1
        truncated(Normal(state[3], 0.1), lower=0.), # sigma2
        truncated(Normal(state[4], 0.1), lower=0.), # sigma3
        MvNormal(state[5], Diagonal([state[2], state[3], state[4]])) # mu
    )
end

AdvancedPS.transition(m::PF, state) = transition_step(m, state)

function obs_density(m::PF, state, t)
    lam_lapse, _, _,_, mu = state
    lapse = logistic(lam_lapse)
    prob = (1 - lapse) * logistic(m.theta.X[t, :]'mu) + lapse * 0.5
    return Bernoulli(prob)
end  

function AdvancedPS.observation(m::PF ,state, t)
    return logpdf(obs_density(m, state, t), y[t])
end

AdvancedPS.isdone(m::PF, t) = t > m.theta.n_trials

n_particles = 20
n_samples = 200
rng = MersenneTwister(2342)

theta0 = Parameters(
    (-9, zeros(3), zeros(3), n_trials, n_cols)
    )

model = PF(theta0)
pgas = AdvancedPS.PGAS(n_particles)
chains = sample(rng, model, pgas, n_samples; progress=true);

There is an issue with the way we handle inner state container in PG and PGAS. For a model like this:

mutable struct Model 
   States::Vector{Float64}
end

PGAS would require the constructor of Model to create an empty array and the sampler would push new states / values into it.
For PG, that's the opposite. We need to start with a pre-allocated array (like zeros(T, N)).

Do we need to target julia 1.3 since the new long term stable version is 1.6 ?

SMC.jl has a SMC implementation with tempering. They claim it works well on problems I am interested in. Would it be possible to implement this with AdvancedPS.jl?

The PG sampler can be improved by adding the ancestor resampling step. It is challenging to implement this for the general case, but perhaps we can consider adding it a a option, e.g.

s = PG(10, 20, :s, as = true)

Here when as equals true, we support only one variable (e.g. :s) to be sampled by PG.

References:

• http://jmlr.org/papers/volume15/lindsten14a/lindsten14a.pdf
• https://papers.nips.cc/paper/5968-particle-gibbs-for-infinite-hidden-markov-models.pdf

I looked into the newly released Threads.@threading construct in Julia 0.5. It seems that adding parallelization feature to SMC is simple - could be done in a few lines of code (see the parallelsmc branch).

However, the threading feature in Julia is still quite fragile (see e.g. here). I've also filed a bug in the Julia repo.

Maybe we should wait until the Julia team fix these threading bugs and revisit to this feature in a few months time.

I had a short conversation with Fredrik Lindsten about SMC algorithms that might be interesting to implement. Apart from PGAS, his suggestions were:

• Our divide and conquer SMC is something that I think would be very well suited for PPL. I have mentioned this to a lot of people (including Lawrence, David Broman, ...) but no one has taken the bait yet (https://arxiv.org/abs/1406.4993).
• Some type of "Twisted SMC" would also be very beneficial for PPL, but that's more of a research topic to figure out how to construct the sampler itself so that it can be implemented in a generic way.
• Last but not least, Christian's VSMC (https://arxiv.org/abs/1705.11140).

This package implements a splittable random number generator based on the standard Random library. It might be helpful to compare against our internal splittable rng.

As discussed in TuringLang/Turing.jl#426, it would be good if Turing would be able to store statistics on the particle weights for all particle-filtering based methods. These statistics could be stored in the ParticleContainer and exposed to the user through new fields in the Chain type.

I have an implementation of Annealed Importance Sampling (AIS) at https://github.com/treigerm/AnnealedIS.jl and @ParadaCarleton kindly suggested that I could merge it into an official Turing package and suggested this one. I was wondering what the best process is to merge my package into this one and if this is the best place.

Skimming through this repository it seems to me that the requirement for AIS and SMC are actually fairly different. AIS technically also has a set of particles but crucially there is no resampling step at observe statements. For a given Turing model AIS requires access to the following functions:

• The log joint density function of the parameters
• The log prior density function of the parameters
  If one wants to use AIS with an HMC transition kernel then we also requires access to the gradient of these functions. Additionally, ideally we want to be able to evaluate these function in parallel on a set of particles. You can see how I currently extract those functions from a Turing model here.

In summary, I have two questions:

• How should I add my implementation of AIS to this package? I.e. should I add it to the contrib folder?
• What is the "best practice" to extract the two functions mentioned above? Am I doing it correctly in my implementation?

Implementation:

• Implement PG and SMC
• Implement PGAS
• Write interface for package with additional functionalities
• Improve exception handling

Testing:

• Test PG and SMC
• Test PGAS
• Test PG and SMC in combination with Turing branch ParticleGibbsExtension
• Test PGAS in combination with Turing branch ParticleGibbsExtension

The README uses the name "AdvancedParticleSamplers.jl" (so maybe that was the initial name of the repo?). IMO that's a lot more descriptive than AdvancedPS.jl, at least I wasn't familiar with the abbreviation PS (in contrast to MH or HMC with their respective packages AdvancedMH and AdvancedHMC). So maybe it could be helpful (also to make it easier to discover this package) to rename it to AdvancedParticleSamplers.jl or something like AdvancedParticleFilters.jl?

Long names shouldn't be a problem for users or developers since one could always define const AdvancedPF = AdvancedParticleFilters at the top of a module or script, if needed.

Differential evolution MCMC is a particle sampler that works well for ABC and other situations in which AD is not feasible. Although it is less efficient than NUTS, it performs better with correlated/irregular posterior distributions compared to many other alternatives. I think it would be a useful addition to Turing. I have a naive implementation in DifferentialEvolutionMCMC that partially follows the AbstractMCMC template. Nonetheless, I think it would benefit from the expertise of Turing developers, who could integrate it better with the Turing ecosystem, make it composable, and optimize performance. If anyone is interested in this algorithm, feel free to reuse any code in DifferentialEvolutionMCMC.

References:
Approximate Bayesian computation with differential evolution