JuliaPOMDP / ParticleFilters.jl

Simple particle filter implementation in Julia - works with POMDPs.jl models or others.
Other
45 stars 16 forks source link

How to use resample? #45

Closed idontgetoutmuch closed 3 years ago

idontgetoutmuch commented 3 years ago

I am trying to use particle filtering on an agent based model.

m = ParticleFilterModel{AgentBasedModel{Union{Grass, Sheep, Wolf},GridSpace{LightGraphs.SimpleGraphs.SimpleGraph{Int64},2,Int64},Agents.var"#by_union#9"{Bool,Bool},Nothing}}(f, g);

# We've already run the model to create some data
zs = convert(Array{Array{Float64, 1}, 1}, [[a, b] for (a, b) in zip(results.count_wolves, results.count_sheep)])
vs = [[a, a] for a in zeros(Float64, 501, 1)]

n = 5
fil = SIRParticleFilter(m, n)

d = Bernoulli(0.5)

Random.seed!(42)

# construct initial belief
b0 = ParticleCollection([
    initialize_model(
        n_sheep = 100 + sum(rand(d,10)),
        n_wolves = 20 + sum(rand(d,10)),
        dims = (25, 25),
        regrowth_time = 30,
        Δenergy_sheep   = exp(log(5)    + 1.0 * randn()),
        Δenergy_wolf    = exp(log(20)   + 1.0 * randn()),
        sheep_reproduce = exp(log(0.2) + 0.1 * randn()),
        wolf_reproduce  = exp(log(0.08) + 0.1 * randn())
    ) for i in 1:n])

df0 = init_model_dataframe(particle(b0,1), mdata)
for i in 1:n collect_model_data!(df0,particle(b0,i),mdata) end

pm = predict(m, b0, vs[1], MersenneTwister(42))

df1 = init_model_dataframe(pm[1], mdata)
for i in 1:n collect_model_data!(df1,pm[i],mdata) end

wm = reweight(m, b0, vs[1], pm, zs[1])

I get weights of

julia> wm = reweight(m, b0, vs[1], pm, zs[1])
5-element Array{Float64,1}:
 0.00026809439279891365
 1.1022573863206847e-18
 4.0029031060636934e-14
 0.00014350058766747782
 2.0188548250289625e-20

so when I resample I should see about 3/4 particles sampled from the particle with weight 0.00026809439279891365 and about 1/2 particles sampled from the particle with weight 0.00014350058766747782.

I can see resample(resampler, bp, rng) but

  1. Where do I get a resampler from and
  2. How do I convert an array to a WeightedParticleBelief

Apologies if this is obvious but I am very new to julia - thanks in advance :-)

idontgetoutmuch commented 3 years ago

This seems to get me some of the way

WeightedParticleBelief(pm,wm,sum(wm))
idontgetoutmuch commented 3 years ago

b1 = resample(ImportanceResampler(n), WeightedParticleBelief(pm,wm,sum(wm)), MersenneTwister(42)) seems to work.