Bayesian Markov-Chain-Monte-Carlo and Claims Data

JuliaActuary is an ecosystem of packages that makes Julia the easiest language to get started for actuarial workflows.
mortalitytables
exposures
experience-analysis
dataframes
tutorial
statistics
bayesian
using MortalityTables
using Turing
using UUIDs
using DataFramesMeta
using MCMCChains, Plots, StatsPlots
using LinearAlgebra
using Pipe
using StatsFuns

Generating fake data

The problem of interest is to look at mortality rates, which are given in terms of exposures (whether or not a life experienced a death in a given year).

We’ll grab some example rates from an insurance table, which has a “selection” component: When someone enters observation, say at age 50, their mortality is path dependent (so for someone who started being observed at 50 will have a different risk/mortality rate at age 55 than someone who started being observed at 45).

Addtionally, there may be additional groups of interest, such as: - high/medium/low risk classification - sex - group (e.g. company, data source, etc.) - type of insurance product offered

The example data will start with only the risk classification above

src = MortalityTables.table("2001 VBT Residual Standard Select and Ultimate - Male Nonsmoker, ANB")
MortalityTable (Insured Lives Mortality):
   Name:
       2001 VBT Residual Standard Select and Ultimate - Male Nonsmoker, ANB
   Fields: 
       (:select, :ultimate, :metadata)
   Provider:
       Society of Actuaries
   mort.SOA.org ID:
       1118
   mort.SOA.org link:
       https://mort.soa.org/ViewTable.aspx?&TableIdentity=1118
   Description:
       2001 Valuation Basic Table (VBT) Residual Standard Select and Ultimate Table -  Male Nonsmoker. Basis: Age Nearest Birthday. Minimum Select Age: 0. Maximum Select Age: 99. Minimum Ultimate Age: 25. Maximum Ultimate Age: 120
n = 10_000

function generate_data_individual(tbl, issue_age=rand(50:55), inforce_years=rand(1:30), risklevel=rand(1:3))
    # risk_factors will scale the "true" parameter up or down
    # we observe the assigned risklevel, but not risk_factor
    risk_factors = [0.7, 1.0, 1.5]
    rf = risk_factors[risklevel]
    deaths = rand(inforce_years) .< (tbl.select[issue_age][issue_age.+inforce_years.-1] .* rf)

    endpoint = if sum(deaths) == 0
        last(inforce_years)
    else
        findfirst(deaths)
    end
    id = uuid1()
    map(1:endpoint) do i
        (
            issue_age=issue_age,
            risklevel=risklevel,
            att_age=issue_age + i - 1,
            death=deaths[i],
            id=id,
        )
    end

end

exposures = vcat([generate_data_individual(src) for _ in 1:n]...) |> DataFrame
108913×5 DataFrame
108888 rows omitted
Row issue_age risklevel att_age death id
Int64 Int64 Int64 Bool UUID
1 54 2 54 false 19d9bd90-2c56-11ef-221d-ab3d14a6f167
2 54 2 55 false 19d9bd90-2c56-11ef-221d-ab3d14a6f167
3 55 1 55 false 19dd329a-2c56-11ef-246b-334b05232974
4 55 1 56 false 19dd329a-2c56-11ef-246b-334b05232974
5 55 1 57 true 19dd329a-2c56-11ef-246b-334b05232974
6 55 2 55 false 19dd32f4-2c56-11ef-1da4-c926e1b3d176
7 55 2 56 false 19dd32f4-2c56-11ef-1da4-c926e1b3d176
8 55 2 57 false 19dd32f4-2c56-11ef-1da4-c926e1b3d176
9 55 2 58 false 19dd32f4-2c56-11ef-1da4-c926e1b3d176
10 55 2 59 false 19dd32f4-2c56-11ef-1da4-c926e1b3d176
11 55 2 60 false 19dd32f4-2c56-11ef-1da4-c926e1b3d176
12 55 2 61 false 19dd32f4-2c56-11ef-1da4-c926e1b3d176
13 55 2 62 false 19dd32f4-2c56-11ef-1da4-c926e1b3d176
⋮ ⋮ ⋮ ⋮ ⋮ ⋮
108902 52 3 56 false 19deb6c4-2c56-11ef-2c57-9b6e03cb9296
108903 52 3 57 false 19deb6c4-2c56-11ef-2c57-9b6e03cb9296
108904 52 3 58 false 19deb6c4-2c56-11ef-2c57-9b6e03cb9296
108905 52 3 59 false 19deb6c4-2c56-11ef-2c57-9b6e03cb9296
108906 52 3 60 false 19deb6c4-2c56-11ef-2c57-9b6e03cb9296
108907 52 3 61 false 19deb6c4-2c56-11ef-2c57-9b6e03cb9296
108908 52 3 62 false 19deb6c4-2c56-11ef-2c57-9b6e03cb9296
108909 52 3 63 false 19deb6c4-2c56-11ef-2c57-9b6e03cb9296
108910 52 3 64 false 19deb6c4-2c56-11ef-2c57-9b6e03cb9296
108911 52 3 65 false 19deb6c4-2c56-11ef-2c57-9b6e03cb9296
108912 52 3 66 false 19deb6c4-2c56-11ef-2c57-9b6e03cb9296
108913 52 3 67 true 19deb6c4-2c56-11ef-2c57-9b6e03cb9296

Two groupings, one with and without the risk level:

data = combine(groupby(exposures, [:issue_age, :att_age])) do subdf
    (exposures=nrow(subdf),
        deaths=sum(subdf.death),
        fraction=sum(subdf.death) / nrow(subdf))
end
180×5 DataFrame
155 rows omitted
Row issue_age att_age exposures deaths fraction
Int64 Int64 Int64 Int64 Float64
1 50 50 1672 35 0.020933
2 50 51 1576 22 0.0139594
3 50 52 1509 28 0.0185553
4 50 53 1412 19 0.0134561
5 50 54 1334 32 0.023988
6 50 55 1257 38 0.0302307
7 50 56 1159 23 0.0198447
8 50 57 1080 22 0.0203704
9 50 58 996 27 0.0271084
10 50 59 921 21 0.0228013
11 50 60 861 27 0.0313589
12 50 61 781 24 0.0307298
13 50 62 709 11 0.0155148
⋮ ⋮ ⋮ ⋮ ⋮ ⋮
169 55 73 220 8 0.0363636
170 55 74 182 12 0.0659341
171 55 75 147 7 0.047619
172 55 76 122 3 0.0245902
173 55 77 102 6 0.0588235
174 55 78 76 5 0.0657895
175 55 79 56 7 0.125
176 55 80 37 5 0.135135
177 55 81 26 1 0.0384615
178 55 82 21 3 0.142857
179 55 83 12 0 0.0
180 55 84 5 0 0.0 
# ╔═╡ 45237199-f8e8-4f61-b644-89ab37c31a5d
data2 = combine(groupby(exposures, [:issue_age, :att_age, :risklevel])) do subdf
    (exposures=nrow(subdf),
        deaths=sum(subdf.death),
        fraction=sum(subdf.death) / nrow(subdf))
end
531×6 DataFrame
506 rows omitted
Row issue_age att_age risklevel exposures deaths fraction
Int64 Int64 Int64 Int64 Int64 Float64
1 50 50 1 550 16 0.0290909
2 50 50 2 601 9 0.014975
3 50 50 3 521 10 0.0191939
4 50 51 1 518 6 0.011583
5 50 51 2 568 7 0.0123239
6 50 51 3 490 9 0.0183673
7 50 52 1 498 6 0.0120482
8 50 52 2 542 12 0.0221402
9 50 52 3 469 10 0.021322
10 50 53 1 474 1 0.0021097
11 50 53 2 507 6 0.0118343
12 50 53 3 431 12 0.0278422
13 50 54 1 450 7 0.0155556
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
520 55 80 2 8 2 0.25
521 55 80 3 4 0 0.0
522 55 81 1 19 0 0.0
523 55 81 2 3 1 0.333333
524 55 81 3 4 0 0.0
525 55 82 1 15 2 0.133333
526 55 82 2 2 0 0.0
527 55 82 3 4 1 0.25
528 55 83 1 9 0 0.0
529 55 83 2 1 0 0.0
530 55 83 3 2 0 0.0
531 55 84 1 5 0 0.0

1: A single binomial parameter model

Estiamte \(p\), the average mortality rate, not accounting for any variation within the population/sample:

@model function mortality(data, deaths)
    p ~ Beta(1, 1)
    for i = 1:nrow(data)
        deaths[i] ~ Binomial(data.exposures[i], p)
    end
end


m1 = mortality(data, data.deaths)

num_chains = 4
4

Sampling from the posterior

We use a No-U-Turn-Sampler (NUTS) technique to sample multile chains at once:

chain = sample(m1, NUTS(), 1000)

plot(chain)
┌ Info: Found initial step size
└   ϵ = 0.0125
Sampling:   4%|█▌                                       |  ETA: 0:00:03Sampling: 100%|█████████████████████████████████████████| Time: 0:00:00

Plotting samples from the posterior

We can see that the sampling of possible posterior parameters doesn’t really fit the data very well since our model was so simplified. The lines represent the posterior binomial probability.

This is saying that for the observed data, if there really is just a single probability p that governs the true process that came up with the data, there’s a pretty narrow range of values it could possibly be:

let
    data_weight = data.exposures ./ sum(data.exposures)
    data_weight = .√(data_weight ./ maximum(data_weight) .* 20)

    p = scatter(
        data.att_age,
        data.fraction,
        markersize=data_weight,
        alpha=0.5,
        label="Experience data point (size indicates relative exposure quantity)",
        xlabel="age",
        ylim=(0.0, 0.25),
        ylabel="mortality rate",
        title="Parametric Bayseian Mortality"
    )

    # show n samples from the posterior plotted on the graph
    n = 300
    ages = sort!(unique(data.att_age))

    for i in 1:n
        p_posterior = sample(chain, 1)[:p][1]
        hline!([p_posterior], label="", alpha=0.1)
    end
    p

end

The posterior mean of p is of course very close to the simple proportoin of claims to exposures:

let
    a = mean(chain, :p)
    b = sum(data.deaths) / sum(data.exposures)
    a, b
end
(0.03054276093985464, 0.030538135943367642)

2. Parametric model

In this example, we utilize a MakehamBeard parameterization because it’s already very similar in form to a logistic function. This is important because our desired output is a probability (ie the probablity of a death at a given age), so the value must be constrained to be in the interval between zero and one.

The prior values for a,b,c, and k are chosen to constrain the hazard (mortality) rate to be between zero and one.

This isn’t an ideal parameterization (e.g. we aren’t including information about the select underwriting period), but is an example of utilizing Bayesian techniques on life experience data.

@model function mortality2(data, deaths)
    a ~ Exponential(0.1)
    b ~ Exponential(0.1)
    c = 0.0
    k ~ truncated(Exponential(1), 1, Inf)

    # use the variables to create a parametric mortality model
    m = MortalityTables.MakehamBeard(; a, b, c, k)

    # loop through the rows of the dataframe to let Turing observe the data 
    # and how consistent the parameters are with the data
    for i = 1:nrow(data)
        age = data.att_age[i]
        q = MortalityTables.hazard(m, age)
        deaths[i] ~ Binomial(data.exposures[i], q)
    end
end
mortality2 (generic function with 2 methods)

Sampling from the posterior

We combine the model with the data and a use use a No-U-Turn-Sampler (NUTS) technique to sample:

m2 = mortality2(data, data.deaths)
chain2 = sample(m2, NUTS(), 1000)
summarize(chain2)
┌ Info: Found initial step size
└   ϵ = 0.0125
Sampling:   2%|▊                                        |  ETA: 0:00:06Sampling:   3%|█▎                                       |  ETA: 0:00:12Sampling:   4%|█▊                                       |  ETA: 0:00:14Sampling:   5%|█▉                                       |  ETA: 0:00:15Sampling:   6%|██▎                                      |  ETA: 0:00:15Sampling:   6%|██▋                                      |  ETA: 0:00:14Sampling:   7%|███                                      |  ETA: 0:00:14Sampling:   8%|███▎                                     |  ETA: 0:00:14Sampling:   9%|███▉                                     |  ETA: 0:00:14Sampling:  11%|████▍                                    |  ETA: 0:00:13Sampling:  12%|█████                                    |  ETA: 0:00:12Sampling:  14%|█████▌                                   |  ETA: 0:00:11Sampling:  16%|██████▌                                  |  ETA: 0:00:10Sampling:  18%|███████▎                                 |  ETA: 0:00:09Sampling:  19%|███████▉                                 |  ETA: 0:00:09Sampling:  21%|████████▊                                |  ETA: 0:00:08Sampling:  23%|█████████▌                               |  ETA: 0:00:08Sampling:  26%|██████████▌                              |  ETA: 0:00:07Sampling:  28%|███████████▎                             |  ETA: 0:00:07Sampling:  30%|████████████▍                            |  ETA: 0:00:06Sampling:  33%|█████████████▍                           |  ETA: 0:00:06Sampling:  34%|██████████████▏                          |  ETA: 0:00:06Sampling:  36%|██████████████▉                          |  ETA: 0:00:05Sampling:  39%|████████████████▏                        |  ETA: 0:00:05Sampling:  42%|█████████████████                        |  ETA: 0:00:05Sampling:  44%|██████████████████▏                      |  ETA: 0:00:04Sampling:  47%|███████████████████▏                     |  ETA: 0:00:04Sampling:  48%|███████████████████▉                     |  ETA: 0:00:04Sampling:  50%|████████████████████▋                    |  ETA: 0:00:04Sampling:  53%|█████████████████████▋                   |  ETA: 0:00:03Sampling:  55%|██████████████████████▍                  |  ETA: 0:00:03Sampling:  57%|███████████████████████▍                 |  ETA: 0:00:03Sampling:  59%|████████████████████████▏                |  ETA: 0:00:03Sampling:  61%|████████████████████████▉                |  ETA: 0:00:03Sampling:  63%|█████████████████████████▉               |  ETA: 0:00:03Sampling:  65%|██████████████████████████▋              |  ETA: 0:00:02Sampling:  67%|███████████████████████████▌             |  ETA: 0:00:02Sampling:  69%|████████████████████████████▎            |  ETA: 0:00:02Sampling:  71%|█████████████████████████████▎           |  ETA: 0:00:02Sampling:  74%|██████████████████████████████▎          |  ETA: 0:00:02Sampling:  76%|███████████████████████████████          |  ETA: 0:00:02Sampling:  77%|███████████████████████████████▊         |  ETA: 0:00:02Sampling:  79%|████████████████████████████████▌        |  ETA: 0:00:01Sampling:  82%|█████████████████████████████████▌       |  ETA: 0:00:01Sampling:  84%|██████████████████████████████████▎      |  ETA: 0:00:01Sampling:  85%|███████████████████████████████████      |  ETA: 0:00:01Sampling:  88%|████████████████████████████████████     |  ETA: 0:00:01Sampling:  90%|████████████████████████████████████▉    |  ETA: 0:00:01Sampling:  92%|█████████████████████████████████████▋   |  ETA: 0:00:01Sampling:  94%|██████████████████████████████████████▋  |  ETA: 0:00:00Sampling:  96%|███████████████████████████████████████▍ |  ETA: 0:00:00Sampling:  98%|████████████████████████████████████████▏|  ETA: 0:00:00Sampling:  99%|████████████████████████████████████████▉|  ETA: 0:00:00Sampling: 100%|█████████████████████████████████████████| Time: 0:00:06
  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯
           a    0.0036    0.0006    0.0000   240.9349   198.5353    0.9998     ⋯
           b    0.0366    0.0031    0.0002   215.0610   171.8594    1.0001     ⋯
           k    2.0919    1.1672    0.0756   269.2057   214.9856    1.0071     ⋯
                                                                1 column omitted

plot(chain2)

Plotting samples from the posterior

We can see that the sampling of possible posterior parameters fits the data well:

let
    data_weight = data.exposures ./ sum(data.exposures)
    data_weight = .√(data_weight ./ maximum(data_weight) .* 20)

    p = scatter(
        data.att_age,
        data.fraction,
        markersize=data_weight,
        alpha=0.5,
        label="Experience data point (size indicates relative exposure quantity)",
        xlabel="age",
        ylim=(0.0, 0.25),
        ylabel="mortality rate",
        title="Parametric Bayseian Mortality"
    )


    # show n samples from the posterior plotted on the graph
    n = 300
    ages = sort!(unique(data.att_age))

    for i in 1:n
        s = sample(chain2, 1)
        a = only(s[:a])
        b = only(s[:b])
        k = only(s[:k])
        c = 0
        m = MortalityTables.MakehamBeard(; a, b, c, k)
        plot!(ages, age -> MortalityTables.hazard(m, age), alpha=0.1, label="")
    end
    p
end

3. Parametric model

This model extends the prior to create a multi-level model. Each risk class (risklevel) gets its own \(a\) paramater in the MakhamBeard model. The prior for \(a_i\) is determined by the hyperparameter \(\bar{a}\).

@model function mortality3(data, deaths)
    risk_levels = length(levels(data.risklevel))
    b ~ Exponential(0.1)
    ā ~ Exponential(0.1)
    a ~ filldist(Exponential(ā), risk_levels)
    c = 0
    k ~ truncated(Exponential(1), 1, Inf)

    # use the variables to create a parametric mortality model

    # loop through the rows of the dataframe to let Turing observe the data 
    # and how consistent the parameters are with the data
    for i = 1:nrow(data)
        risk = data.risklevel[i]

        m = MortalityTables.MakehamBeard(; a=a[risk], b, c, k)
        age = data.att_age[i]
        q = MortalityTables.hazard(m, age)
        deaths[i] ~ Binomial(data.exposures[i], q)
    end
end
mortality3 (generic function with 2 methods)

Instantiate model with the data, sample, and summarize:

m3 = mortality3(data2, data2.deaths)
chain3 = sample(m3, NUTS(), 1000)
summarize(chain3)
┌ Info: Found initial step size
└   ϵ = 0.0125
Sampling:   1%|▋                                        |  ETA: 0:00:17Sampling:   2%|▊                                        |  ETA: 0:00:19Sampling:   3%|█▏                                       |  ETA: 0:00:34Sampling:   3%|█▎                                       |  ETA: 0:00:52Sampling:   4%|█▌                                       |  ETA: 0:01:06Sampling:   4%|█▊                                       |  ETA: 0:01:03Sampling:   5%|█▉                                       |  ETA: 0:01:06Sampling:   5%|██▏                                      |  ETA: 0:01:08Sampling:   6%|██▎                                      |  ETA: 0:01:10Sampling:   6%|██▌                                      |  ETA: 0:01:09Sampling:   6%|██▋                                      |  ETA: 0:01:10Sampling:   7%|██▉                                      |  ETA: 0:01:13Sampling:   7%|███                                      |  ETA: 0:01:11Sampling:   8%|███▎                                     |  ETA: 0:01:09Sampling:   8%|███▌                                     |  ETA: 0:01:08Sampling:   9%|███▉                                     |  ETA: 0:01:02Sampling:  10%|████                                     |  ETA: 0:01:00Sampling:  11%|████▍                                    |  ETA: 0:00:56Sampling:  11%|████▋                                    |  ETA: 0:00:54Sampling:  12%|████▊                                    |  ETA: 0:00:53Sampling:  13%|█████▏                                   |  ETA: 0:00:50Sampling:  14%|█████▌                                   |  ETA: 0:00:47Sampling:  14%|█████▉                                   |  ETA: 0:00:45Sampling:  15%|██████▍                                  |  ETA: 0:00:42Sampling:  16%|██████▋                                  |  ETA: 0:00:40Sampling:  17%|██████▉                                  |  ETA: 0:00:40Sampling:  17%|███████                                  |  ETA: 0:00:39Sampling:  18%|███████▎                                 |  ETA: 0:00:39Sampling:  18%|███████▌                                 |  ETA: 0:00:38Sampling:  19%|███████▋                                 |  ETA: 0:00:38Sampling:  19%|███████▉                                 |  ETA: 0:00:37Sampling:  20%|████████▎                                |  ETA: 0:00:36Sampling:  20%|████████▍                                |  ETA: 0:00:36Sampling:  21%|████████▊                                |  ETA: 0:00:34Sampling:  22%|█████████                                |  ETA: 0:00:34Sampling:  22%|█████████▏                               |  ETA: 0:00:33Sampling:  23%|█████████▌                               |  ETA: 0:00:32Sampling:  24%|█████████▊                               |  ETA: 0:00:32Sampling:  24%|█████████▉                               |  ETA: 0:00:32Sampling:  25%|██████████▍                              |  ETA: 0:00:30Sampling:  26%|██████████▌                              |  ETA: 0:00:30Sampling:  26%|██████████▊                              |  ETA: 0:00:30Sampling:  27%|███████████▏                             |  ETA: 0:00:29Sampling:  28%|███████████▌                             |  ETA: 0:00:28Sampling:  28%|███████████▋                             |  ETA: 0:00:27Sampling:  29%|████████████                             |  ETA: 0:00:26Sampling:  30%|████████████▎                            |  ETA: 0:00:26Sampling:  31%|████████████▋                            |  ETA: 0:00:26Sampling:  31%|████████████▊                            |  ETA: 0:00:25Sampling:  32%|█████████████▎                           |  ETA: 0:00:25Sampling:  33%|█████████████▋                           |  ETA: 0:00:24Sampling:  34%|██████████████                           |  ETA: 0:00:24Sampling:  35%|██████████████▍                          |  ETA: 0:00:23Sampling:  35%|██████████████▌                          |  ETA: 0:00:23Sampling:  36%|██████████████▉                          |  ETA: 0:00:22Sampling:  37%|███████████████▏                         |  ETA: 0:00:22Sampling:  37%|███████████████▎                         |  ETA: 0:00:22Sampling:  38%|███████████████▌                         |  ETA: 0:00:21Sampling:  39%|███████████████▉                         |  ETA: 0:00:21Sampling:  39%|████████████████▏                        |  ETA: 0:00:21Sampling:  40%|████████████████▌                        |  ETA: 0:00:20Sampling:  41%|████████████████▊                        |  ETA: 0:00:20Sampling:  42%|█████████████████▎                       |  ETA: 0:00:19Sampling:  43%|█████████████████▋                       |  ETA: 0:00:19Sampling:  44%|██████████████████                       |  ETA: 0:00:18Sampling:  45%|██████████████████▍                      |  ETA: 0:00:18Sampling:  46%|██████████████████▊                      |  ETA: 0:00:17Sampling:  47%|███████████████████▏                     |  ETA: 0:00:17Sampling:  48%|███████████████████▌                     |  ETA: 0:00:16Sampling:  49%|████████████████████▎                    |  ETA: 0:00:15Sampling:  52%|█████████████████████▎                   |  ETA: 0:00:14Sampling:  54%|██████████████████████▏                  |  ETA: 0:00:13Sampling:  55%|██████████████████████▌                  |  ETA: 0:00:13Sampling:  56%|███████████████████████                  |  ETA: 0:00:13Sampling:  57%|███████████████████████▍                 |  ETA: 0:00:12Sampling:  57%|███████████████████████▌                 |  ETA: 0:00:12Sampling:  58%|███████████████████████▉                 |  ETA: 0:00:12Sampling:  59%|████████████████████████▏                |  ETA: 0:00:12Sampling:  59%|████████████████████████▎                |  ETA: 0:00:12Sampling:  60%|████████████████████████▌                |  ETA: 0:00:11Sampling:  60%|████████████████████████▋                |  ETA: 0:00:11Sampling:  61%|█████████████████████████                |  ETA: 0:00:11Sampling:  62%|█████████████████████████▍               |  ETA: 0:00:11Sampling:  62%|█████████████████████████▋               |  ETA: 0:00:10Sampling:  63%|██████████████████████████               |  ETA: 0:00:10Sampling:  64%|██████████████████████████▎              |  ETA: 0:00:10Sampling:  65%|██████████████████████████▋              |  ETA: 0:00:10Sampling:  66%|███████████████████████████              |  ETA: 0:00:09Sampling:  66%|███████████████████████████▏             |  ETA: 0:00:09Sampling:  67%|███████████████████████████▍             |  ETA: 0:00:09Sampling:  67%|███████████████████████████▌             |  ETA: 0:00:09Sampling:  68%|███████████████████████████▉             |  ETA: 0:00:09Sampling:  69%|████████████████████████████▏            |  ETA: 0:00:09Sampling:  70%|████████████████████████████▌            |  ETA: 0:00:08Sampling:  70%|████████████████████████████▉            |  ETA: 0:00:08Sampling:  71%|█████████████████████████████▏           |  ETA: 0:00:08Sampling:  72%|█████████████████████████████▌           |  ETA: 0:00:08Sampling:  72%|█████████████████████████████▋           |  ETA: 0:00:08Sampling:  73%|██████████████████████████████           |  ETA: 0:00:07Sampling:  74%|██████████████████████████████▎          |  ETA: 0:00:07Sampling:  75%|██████████████████████████████▋          |  ETA: 0:00:07Sampling:  76%|███████████████████████████████          |  ETA: 0:00:07Sampling:  76%|███████████████████████████████▏         |  ETA: 0:00:06Sampling:  77%|███████████████████████████████▋         |  ETA: 0:00:06Sampling:  77%|███████████████████████████████▊         |  ETA: 0:00:06Sampling:  78%|████████████████████████████████         |  ETA: 0:00:06Sampling:  79%|████████████████████████████████▎        |  ETA: 0:00:06Sampling:  80%|████████████████████████████████▊        |  ETA: 0:00:05Sampling:  81%|█████████████████████████████████▏       |  ETA: 0:00:05Sampling:  81%|█████████████████████████████████▎       |  ETA: 0:00:05Sampling:  82%|█████████████████████████████████▋       |  ETA: 0:00:05Sampling:  83%|██████████████████████████████████       |  ETA: 0:00:04Sampling:  84%|██████████████████████████████████▎      |  ETA: 0:00:04Sampling:  84%|██████████████████████████████████▋      |  ETA: 0:00:04Sampling:  85%|███████████████████████████████████      |  ETA: 0:00:04Sampling:  86%|███████████████████████████████████▏     |  ETA: 0:00:04Sampling:  87%|███████████████████████████████████▋     |  ETA: 0:00:03Sampling:  88%|████████████████████████████████████     |  ETA: 0:00:03Sampling:  88%|████████████████████████████████████▏    |  ETA: 0:00:03Sampling:  89%|████████████████████████████████████▍    |  ETA: 0:00:03Sampling:  89%|████████████████████████████████████▌    |  ETA: 0:00:03Sampling:  90%|████████████████████████████████████▊    |  ETA: 0:00:03Sampling:  90%|█████████████████████████████████████▏   |  ETA: 0:00:02Sampling:  91%|█████████████████████████████████████▌   |  ETA: 0:00:02Sampling:  92%|█████████████████████████████████████▋   |  ETA: 0:00:02Sampling:  92%|█████████████████████████████████████▉   |  ETA: 0:00:02Sampling:  93%|██████████████████████████████████████▎  |  ETA: 0:00:02Sampling:  94%|██████████████████████████████████████▌  |  ETA: 0:00:02Sampling:  95%|██████████████████████████████████████▉  |  ETA: 0:00:01Sampling:  95%|███████████████████████████████████████  |  ETA: 0:00:01Sampling:  96%|███████████████████████████████████████▎ |  ETA: 0:00:01Sampling:  97%|███████████████████████████████████████▋ |  ETA: 0:00:01Sampling:  98%|████████████████████████████████████████ |  ETA: 0:00:01Sampling:  98%|████████████████████████████████████████▏|  ETA: 0:00:01Sampling:  98%|████████████████████████████████████████▍|  ETA: 0:00:00Sampling:  99%|████████████████████████████████████████▌|  ETA: 0:00:00Sampling:  99%|████████████████████████████████████████▊|  ETA: 0:00:00Sampling: 100%|█████████████████████████████████████████| Time: 0:00:25
  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯
           b    0.0447    0.0053    0.0012    33.2602    21.5493    1.0398     ⋯
           ā    0.0065    0.0071    0.0005    54.2101    85.6676    1.0458     ⋯
        a[1]    0.0016    0.0004    0.0001    34.7822    21.6119    1.0297     ⋯
        a[2]    0.0022    0.0005    0.0001    34.5119    18.3744    1.0348     ⋯
        a[3]    0.0035    0.0008    0.0001    37.7816    42.2177    1.0315     ⋯
           k    2.8153    1.9264    0.4387    38.2738    20.6857    1.0268     ⋯
                                                                1 column omitted

let data = data2

    data_weight = data.exposures ./ sum(data.exposures)
    data_weight = .√(data_weight ./ maximum(data_weight) .* 20)
    color_i = data.risklevel

    p = scatter(
        data.att_age,
        data.fraction,
        markersize=data_weight,
        alpha=0.5,
        color=color_i,
        label="Experience data point (size indicates relative exposure quantity)",
        xlabel="age",
        ylim=(0.0, 0.25),
        ylabel="mortality rate",
        title="Parametric Bayseian Mortality"
    )


    # show n samples from the posterior plotted on the graph
    n = 100

    ages = sort!(unique(data.att_age))
    for r in 1:3
        for i in 1:n
            s = sample(chain3, 1)
            a = only(s[Symbol("a[$r]")])
            b = only(s[:b])
            k = only(s[:k])
            c = 0
            m = MortalityTables.MakehamBeard(; a, b, c, k)
            if i == 1
                plot!(ages, age -> MortalityTables.hazard(m, age), label="risk level $r", alpha=0.2, color=r)
            else
                plot!(ages, age -> MortalityTables.hazard(m, age), label="", alpha=0.2, color=r)
            end
        end
    end
    p
end