Approximate Bayesian Computation

& writing performant Python code

Patrick J. Laub and Pierre-Olivier Goffard

https://slides.com/plaub/abc & https://arxiv.org/abs/2007.03833

Motivation

Have a random number of claims \(N \sim p_N( \,\cdot\, ; \boldsymbol{\theta}_{\mathrm{freq}} )\) and the claim sizes \(U_1, \dots, U_N \overset{\mathrm{i.i.d.}}{\sim} f_U( \,\cdot\, ; \boldsymbol{\theta}_{\mathrm{sev}} )\).

We aggregate them somehow, like:

• aggregate claims: \(X = \sum_{i=1}^N U_i \)
• maximum claims: \(X = \max_{i=1}^N U_i \)
• stop-loss: \(X = ( \sum_{i=1}^N U_i - c )_+ \).

Question: Given a sample \(X_1, \dots, X_n\) of the summaries, what is the \(\boldsymbol{\theta} = (\boldsymbol{\theta}_{\mathrm{freq}}, \boldsymbol{\theta}_{\mathrm{sev}})\) which explains them?

Easier question: Given \((X_1, N_1), \dots, (X_n, N_n)\) summaries & counts, what is \(\boldsymbol{\theta}\)?

E.g. a reinsurance contract

Likelihoods

For simple rv's we know their likelihood (normal, exponential, gamma, etc.).

 

When simple rv's are combined, the resulting thing rarely has a likelihood.

$$ X_1, X_2 \overset{\mathrm{i.i.d.}}{\sim} f_X(\,\cdot\,) \Rightarrow X_1 + X_2 \sim ~ \texttt{Unknown Likelihood}! $$

For a sample of \(n\) i.i.d. observations the joint likelihood is

$$ p_{\boldsymbol{X}}(\boldsymbol{x} \mid \boldsymbol{\theta}) = \prod_{i=1}^n p_{X_i}(x_i; \boldsymbol{\theta}) \,. $$

If \(n\) increases, then \(p_{\boldsymbol{X}}(\boldsymbol{x} \mid \boldsymbol{\theta}) = \prod \text{Small things} \overset{\dagger}{=} 0\), or just takes a long time to compute, then \(\texttt{Intractable Likelihood}\)!

Bayesian statistics

Markov chain Monte Carlo

Exact matching algorithm

Given some observations \(\boldsymbol{x}_{\text{obs}}\), repeat:

• generate a potential parameter from the prior distribution \(\boldsymbol{\theta}^{\ast} \sim \pi(\boldsymbol{\theta})\);
• simulate some 'fake data' \(\boldsymbol{x}^{\ast}\) from the model \(\boldsymbol{\theta}^{\ast}\);
• if \( \boldsymbol{x}_{\text{obs}} = \boldsymbol{x}^{\ast}\), then store \(\boldsymbol{\theta}^{\ast}\).

The resulting \(\boldsymbol{\theta}^{\ast}\)s are an i.i.d. sample from the posterior \(\pi(\theta \mid \boldsymbol{x}_{\text{obs}})\).

ABC acceptance–rejection algorithm

Given some observations \(\boldsymbol{x}_{\text{obs}}\), to generate \(R\) samples from the posterior:

• For \(i = 1\), \(\dots\), \(R\):
  • Repeat:
    • generate a potential parameter from the prior distribution \(\boldsymbol{\theta}^{\ast} \sim \pi(\boldsymbol{\theta})\);
    • simulate some 'fake data' \(\boldsymbol{x}^{\ast}\) from the model \(\boldsymbol{\theta}^{\ast}\);
  • Until \( \lVert \boldsymbol{x}_{\text{obs}}-\boldsymbol{x}^{\ast} \rVert \leq \epsilon\)
  • Store \(\boldsymbol{\theta}^{\ast}_i = \boldsymbol{\theta}^{\ast}\)

Inputs: \(\lVert \,\cdot\, \rVert\) denotes a distance measure, and \(\epsilon \ge 0\) is an acceptance threshold.

This samples from the approximate posterior a.k.a. the ABC posterior

$$ \pi_\epsilon(\boldsymbol{\theta} \mid \boldsymbol{x}_{\text{obs}}) \propto \pi(\theta) \times \mathbb{P}\bigl( \lVert \boldsymbol{x}_{\text{obs}}-\boldsymbol{x}^{\ast} \rVert \leq \epsilon \text{ where } \boldsymbol{x}^{\ast} \sim \boldsymbol{\theta} \bigr) . $$

ABC Sequential Monte Carlo

Main idea: Take a bad approximation and make a better one.

Say we have samples from a bad approximation \(\pi_3(\boldsymbol{\theta} \mid \boldsymbol{x}_{\text{obs}})\) and want to improve them to be samples from a better approximation \(\pi_1(\boldsymbol{\theta} \mid \boldsymbol{x}_{\text{obs}})\)?

• Fit a kernel density estimate \(K\) to the samples from the bad distribution \(\pi_3(\boldsymbol{\theta} \mid \boldsymbol{x}_{\text{obs}})\)
• For \(i = 1\), \(\dots\), \(R\):
  • Repeat:
    • generate a potential parameter \(\boldsymbol{\theta}^{\ast}\) from the KDE \(K\);
    • simulate some 'fake data' \(\boldsymbol{x}^{\ast}\) from the model \(\boldsymbol{\theta}^{\ast}\);
  • Until \( \lVert \boldsymbol{x}_{\text{obs}}-\boldsymbol{x}^{\ast} \rVert \leq 1\)
  • Store \(\boldsymbol{\theta}^{\ast}_i = \boldsymbol{\theta}^{\ast}\) along with a weight \(w_i = \frac{ \pi(\boldsymbol{\theta}^{\ast}) }{ K(\boldsymbol{\theta}^{\ast}) }\)

Repeat this to get the ABC Sequential Monte Carlo algorithm, a.k.a. Population Monte Carlo.

Dependent Poisson – Exponential variables

ABC posterior vs true posterior

Observe \(X = \sum_{i=1}^N U_i\) where \(N \sim \mathsf{Poisson}(\lambda=4)\), \(U_i \mid N \sim \mathsf{DepExp}(\beta=2, \delta=0.2)\)

ABC posteriors (KDEs) based on 50 \(X\)'s and on 250 \(X\)'s.

True posteriors (histograms) based on all the \(N\)'s and \(U_i\)'s used to make the 50 \(X\)'s and the 250 \(X\)'s.

Code demonstration

Negative Binomial – Weibull variables

• frequency has negative binomial distribution \(N \sim \mathsf{NegBin}(\alpha=4, p=\frac23)\),
• severities have Weibull distributions \(U_i \sim \mathsf{Weibull}(k=\frac13, \beta=1)\),
• global stop-loss summary \(X = ( \sum_{i=1}^N U_i - 6 )_+\)

ABC posteriors based on 50 (11 non-zero) \(X\)'s and on 250 (69 non-zero) \(X\)'s with uniform priors.

Negative Binomial – Weibull variables: also observe the frequency

Observe \(X\) and \(N\) where \(X = ( \sum_{i=1}^N U_i - 6 )_+\), \(N \sim \mathsf{NegBin}(\alpha=4, p=\frac23)\) and \(U_i \sim \mathsf{Weibull}(k=\frac13, \beta=1)\).

ABC posteriors based on 50 (11 non-zero) \(X\)'s and on 250 (69 non-zero) \(X\)'s with uniform priors.

Same figure zoomed in...

Negative Binomial – Gamma variables (misspecified!)

• frequency has negative binomial distribution \(N \sim \mathsf{NegBin}(\alpha=4, p=\frac23)\),
• severities are fit with \(U_i \sim \mathsf{Gamma}(r, m)\) but generated by \(\mathsf{Weibull}(k=\frac13, \beta=1)\)
• global stop-loss summary \(X = ( \sum_{i=1}^N U_i - 6 )_+\)

ABC posteriors based on 50 (11 non-zero) \(X\)'s and on 250 (69 non-zero) \(X\)'s with uniform priors.

Model selection

Given some observations \(\boldsymbol{x}_{\text{obs}}\), to generate \(R\) samples from the posterior:

• For \(i = 1\), \(\dots\), \(R\):
  • Repeat:
    • choose a model \(m^{\ast} \sim \pi(m)\)
    • generate a potential parameter from the prior distribution of that model \(\boldsymbol{\theta}^{\ast} \sim \pi(\boldsymbol{\theta} \mid m)\);
    • simulate some 'fake data' \(\boldsymbol{x}^{\ast}\) from that parameter \(\boldsymbol{\theta}^{\ast}\);
  • Until \( \lVert \boldsymbol{x}_{\text{obs}}-\boldsymbol{x}^{\ast} \rVert \leq \epsilon\)
  • Store \(m^{\ast}_i = m^{\ast}\)

Then the \(m^{\ast}\)s are samples from \(\pi(m \mid \boldsymbol{x}_{\text{obs}})\).

Model posterior: Weibull (✓) vs gamma (✗)?

Writing faster code

Contributions

• Modest theoretical contribution, looking at convergence of the approximate posterior to the true posterior in the case of mixed data (e.g. truncated data).
• An optimized parallel Python package abcre ('re' for reinsurance), on Github

- What ABC can do, when to use it (as a last resort)

- ABC turns a statistics problem into a programming problem

- Profile your code (snakeviz)

- Vectorise your code (numpy)

- Compile your code (numba)

- Parallelise your code (joblib), don't overdo it