# Poisson changepoint model: UK Coal Mining Disasters

**URL:** <https://discourse.julialang.org/t/poisson-changepoint-model-uk-coal-mining-disasters/57470>\
**Category:** Probabilistic Programming\
**Tags:** question\
**Created:** [March 18, 2021, 12:58pm UTC](https://discourse.julialang.org/t/poisson-changepoint-model-uk-coal-mining-disasters/57470 "2021-03-18T12:58:35Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![gsinha](https://avatars.discourse-cdn.com/v4/letter/g/57b2e6/32.png) [@gsinha](https://discourse.julialang.org/u/gsinha)\
**Post date:** [March 18, 2021, 12:58pm UTC](https://discourse.julialang.org/t/poisson-changepoint-model-uk-coal-mining-disasters/57470/1 "2021-03-18T12:58:36Z")

</div>

Hi:

In an effort to learn Turing, I set out to replicate the [UK Coal Mining disasters example](https://colcarroll.github.io/pymc3/notebooks/getting_started.html). I think I’ve managed to successfully replicate the results - this is what I did.

```julia
using Turing
using Distributions
using Distributed 
using RollingFunctions

using MCMCChains, Plots, StatsPlots

use_interp_data = false
if use_interp_data 
    disasters = [
        4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
        3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
        2, 2, 3, 4, 2, 1, 3, missing, 2, 1, 1, 1, 1, 3, 0, 0,
        1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
        0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
        3, 3, 1, missing, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
        0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1
    ]
else 
    disasters = [
        4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
        3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
        2, 2, 3, 4, 2, 1, 3, 0, 2, 1, 1, 1, 1, 3, 0, 0,
        1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
        0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
        3, 3, 1, 0, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
        0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1
    ]
end

years = collect(1851:1961)

scatter(years, disasters, markertype=:circle, markersize=8, 
    alpha=0.4, xlabel="Year", legend=:false, ylabel="Disaster count")

@model function disasters_model(x, y)
    s ~ Uniform(1851, 1961)
    λ₀ ~ Exponential(1.)
    λ₁ ~ Exponential(1.)

    λ = map(z -> (z <= s) ? λ₀ : λ₁, x)
    y .~ Poisson.(λ)
end

n_samples = 2000
n_chains = 4
delta = .80
n_adapt = 1000
config = NUTS(n_adapt, delta)

addprocs(4)

@time chain = sample(disasters_model(years, disasters), 
    config, MCMCThreads(), n_samples, n_chains)

```

The results are below:

```julia
┌ Info: Found initial step size
└ ϵ = 0.4
┌ Info: Found initial step size
└ ϵ = 0.2
┌ Info: Found initial step size
└ ϵ = 0.2
┌ Info: Found initial step size
└ ϵ = 0.2
129.208340 seconds (816.63 M allocations: 143.230 GiB, 33.76% gc time)
Chains MCMC chain (2000×15×4 Array{Float64,3}):

Iterations = 1:2000
Thinning interval = 1
Chains = 1, 2, 3, 4
Samples per chain = 2000
parameters = s, λ₀, λ₁
internals = acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, is_accept, log_density, lp, max_hamiltonian_energy_error, n_steps, nom_step_size, numerical_error, step_size, tree_depth

Summary Statistics
  parameters mean std naive_se mcse ess rhat 
      Symbol Float64 Float64 Float64 Float64 Float64 Float64 

           s 1903.7427 24.8248 0.2775 2.7897 16.0908 15.9049
          λ₀ 2.7723 0.5560 0.0062 0.0590 19.1846 2.2552
          λ₁ 0.7362 0.2725 0.0030 0.0287 18.6696 2.6390

Quantiles
  parameters 2.5% 25.0% 50.0% 75.0% 97.5% 
      Symbol Float64 Float64 Float64 Float64 Float64 

           s 1886.1786 1889.0095 1889.9111 1908.5791 1947.3571
          λ₀ 1.7485 2.3498 2.8690 3.1567 3.6972
          λ₁ 0.1467 0.6906 0.8224 0.9219 1.0905

```

I have a couple of questions for this forum:

- is the code structure efficient? the PyMC3 example runs in roughly 17 seconds versus the 130 seconds here.
- how does one deal with the missing values in the data; I couldn’t figure out how to do it so ended up replacing them with 0 (I wouldn’t do that in real life but just wanted to get something going).
- any other comments?

Thanks for your help.

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [March 18, 2021, 1:09pm UTC](https://discourse.julialang.org/t/poisson-changepoint-model-uk-coal-mining-disasters/57470/2 "2021-03-18T13:09:15Z")

</div>

Not sure of Turing specifics here, but `map` will allocate at each iteration, leading to lots of overhead. You might try a `MappedArray` here:

```julia

julia> using MappedArrays, BenchmarkTools

julia> x = randn(10000);

julia> @btime map(z -> (z ≤ 0.5 ? 2.0 : 3.0), $x);
  7.106 μs (2 allocations: 78.20 KiB)

julia> @btime mappedarray(z -> (z ≤ 0.5 ? 2.0 : 3.0), $x);
  1.182 ns (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [March 18, 2021, 1:25pm UTC](https://discourse.julialang.org/t/poisson-changepoint-model-uk-coal-mining-disasters/57470/3 "2021-03-18T13:25:32Z")

</div>

Uh, why is there missing data? Is this something added for education purposes, I think I remember there are no missing years in the original data set [https://academic.oup.com/biomet/article/66/1/191/224063](https://academic.oup.com/biomet/article/66/1/191/224063)

---

<div class="post-metadata">

**Author:** ![gsinha](https://avatars.discourse-cdn.com/v4/letter/g/57b2e6/32.png) [@gsinha](https://discourse.julialang.org/u/gsinha)\
**Post date:** [March 18, 2021, 1:41pm UTC](https://discourse.julialang.org/t/poisson-changepoint-model-uk-coal-mining-disasters/57470/4 "2021-03-18T13:41:24Z")

</div>

Good question - I don’t have access to Biometrika so will have to take your word for it. I simply replicated the PyMC3 example dataset for this purpose.

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [March 18, 2021, 2:41pm UTC](https://discourse.julialang.org/t/poisson-changepoint-model-uk-coal-mining-disasters/57470/5 "2021-03-18T14:41:00Z")

</div>

Yes, I see that its about Turing and not about coal mines, I was just puzzled. We had a look at the same data set for [PointProcessInference.jl](https://github.com/mschauer/PointProcessInference.jl#example-2), but there we use event times\* instead of aggregate events.

\*from R:

```nohighlight
> install.packages("boot")
> library(boot)
> coal

```

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [March 18, 2021, 2:44pm UTC](https://discourse.julialang.org/t/poisson-changepoint-model-uk-coal-mining-disasters/57470/6 "2021-03-18T14:44:00Z")

</div>

You can get it from `RDatasets`:

```julia
julia> using RDatasets

julia> coal = dataset("boot","coal")
191×1 DataFrame
 Row │ Date    
     │ Float64 
─────┼─────────
   1 │ 1851.2
   2 │ 1851.63
   3 │ 1851.97
   4 │ 1851.97
   5 │ 1852.31
   6 │ 1852.35
   7 │ 1852.36
...

```

---

<div class="post-metadata">

**Author:** ![gsinha](https://avatars.discourse-cdn.com/v4/letter/g/57b2e6/32.png) [@gsinha](https://discourse.julialang.org/u/gsinha)\
**Post date:** [March 18, 2021, 3:48pm UTC](https://discourse.julialang.org/t/poisson-changepoint-model-uk-coal-mining-disasters/57470/7 "2021-03-18T15:48:48Z")

</div>

Thanks for the suggestion. That brought it down to roughly 85 seconds, substantial progress.

---

<div class="post-metadata">

**Author:** ![BradGroff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bradgroff/32/21156_2.png) [@BradGroff](https://discourse.julialang.org/u/BradGroff)\
**Post date:** [March 18, 2021, 5:54pm UTC](https://discourse.julialang.org/t/poisson-changepoint-model-uk-coal-mining-disasters/57470/8 "2021-03-18T17:54:36Z")

</div>

I didn’t know about the `MappedArray` library, thanks!

So to be explicit, the difference between array and mappedarray in terms of performance here is that on every iteration calling `map` will explicitly create a new array thus requiring memory allocation of that array, and then for lookup each z value one by one when sampling the line `y .~ Poisson.(\lambda)`. So, array allocations = O(n\_samples), function calls = O(n\_samples \* len(y)). I suppose there’s also garbage collection for each of those n\_samples arrays.

In contrast, mapped array skips the allocation step (and hence garbage collection) and just calculates z when each element is accessed. Same number of function calls in each case since, but the array allocation is skipped entirely.

Is there also something else going on here though? Some naive back-of-the-envelop math: you showed you can save a few microseconds allocating an array of size 10k. Accumulated over 2000 samples (ignoring that this is much smaller data), that adds up to only a handful of milliseconds. Even if that isn’t counting gc, that will still be far less than the 45 seconds claimed here. Obviously different machines but still a big gap. Am I missing something?

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [March 18, 2021, 5:57pm UTC](https://discourse.julialang.org/t/poisson-changepoint-model-uk-coal-mining-disasters/57470/9 "2021-03-18T17:57:58Z")

</div>

Right, this is because a `MappedArray` is lazy. The comparison I showed doesn’t include actually computing the elements. So maybe this is better:

```julia
julia> @btime begin
       y = map(z -> (z ≤ 0.5 ? 2.0 : 3.0), $x)
       sum(y)
       end
  5.884 μs (2 allocations: 78.20 KiB)
23126.0

julia> @btime begin
       y = mappedarray(z -> (z ≤ 0.5 ? 2.0 : 3.0), $x)
       sum(y)
       end
  1.534 μs (0 allocations: 0 bytes)
23126.0

```

Also, note that it’s not always a win. If the computed values are accessed multiple times then it’s often better to allocate.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [March 18, 2021, 6:05pm UTC](https://discourse.julialang.org/t/poisson-changepoint-model-uk-coal-mining-disasters/57470/10 "2021-03-18T18:05:14Z")

</div>

> [@gsinha](#):
>
> ```julia
> λ = map(z -> (z <= s) ? λ₀ : λ₁, x)
> y .~ Poisson.(λ)
> 
> ```

Can both those lines be combined in a single `mappedarray` to avoid the allocation alltogether? Does Turing buy `y ~ mappedarray`?

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [March 19, 2021, 8:32pm UTC](https://discourse.julialang.org/t/poisson-changepoint-model-uk-coal-mining-disasters/57470/11 "2021-03-19T20:32:36Z")

</div>

> [@baggepinnen](#):
>
> Can both those lines be combined in a single `mappedarray` to avoid the allocation alltogether?

Good point. Maybe this?

```julia
y .~ mappedarray(z -> (z <= s) ? Poisson(λ₀) : Poisson(λ₁), x)

```