---
title: "Hierarchical Customer Lifetime Value"
author: "Brian Callander"
date: "2019-05-05"
tags: customer lifetime value, recency frequency, hierarchical model, centred parameterisation, non-centred parameterisation, prior-predictive distribution, stan, e-bfmi, energy 
tldr: "We extend the basic CLV model to use hierarchical priors in two different ways: centred and non-centred parameterisations. I'm not aware of any other HMC-based implementations of this hierarchical CLV model, so we'll run some basic tests to check it's doing the right thing. More specifically, we'll fit it to a dataset drawn from the prior predictive distribution. The resulting fits pass the main diagnostic tests and the 90% posterior intervals capture about 91% of the true parameter values."
always_allow_html: yes
output: 
  md_document:
    variant: markdown
    preserve_yaml: yes
---

In a [previous post](./pareto-nbd.html), we described how a model of customer lifetime value (CLV) works, implemented it in Stan, and fit the model to simulated data.  In this post, we'll extend the model to use hierarchical priors in two different ways: [centred and non-centred](https://mc-stan.org/docs/2_18/stan-users-guide/reparameterization-section.html) parameterisations. I'm not aware of any other HMC-based implementations of this hierarchical CLV model, so we'll run some basic tests to check it's doing the right thing. More specifically, we'll fit it to a dataset drawn from the prior predictive distribution. The resulting fits pass the main diagnostic tests and the 90% posterior intervals capture about 91% of the true parameter values. 

<!--more-->

A word of warning: I came across a number of examples where the model showed severe [E-BFMI](https://mc-stan.org/misc/warnings.html#bfmi-low) problems. This seems to happen when the parameters $\mu$ (inverse expected lifetime) and $\lambda$ (expected purchase rate) are fairly similar, but I haven't pinned down a solid reason for these energy problems yet. I suspect this has something to do with the difficulty in distinguishing a short lifetime from a low purchase rate. This is a topic we'll leave for a future post.

We'll work with raw stan, which can be a bit fiddly sometimes. To keep this post within the limits of readability, I'll define some custom functions to simplify the process. Check out the [full code](./recency_frequency.Rmd) to see the details of these functions.

<div>$\DeclareMathOperator{\dbinomial}{Binomial}\DeclareMathOperator{\dbernoulli}{Bernoulli}\DeclareMathOperator{\dpoisson}{Poisson}\DeclareMathOperator{\dnormal}{Normal}\DeclareMathOperator{\dt}{t}\DeclareMathOperator{\dcauchy}{Cauchy}\DeclareMathOperator{\dexponential}{Exp}\DeclareMathOperator{\duniform}{Uniform}\DeclareMathOperator{\dgamma}{Gamma}\DeclareMathOperator{\dinvgamma}{InvGamma}\DeclareMathOperator{\invlogit}{InvLogit}\DeclareMathOperator{\logit}{Logit}\DeclareMathOperator{\ddirichlet}{Dirichlet}\DeclareMathOperator{\dbeta}{Beta}$</div>

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  cache = TRUE,
  # dev = "svglite",
  echo = TRUE,
  comment = NA,
  message = FALSE,
  warning = FALSE,
  error = TRUE
)

library(tidyverse)
library(scales)

library(kableExtra)
library(here)

library(rstan)
library(tidybayes)
library(bayesplot)


rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

theme_set(theme_bw())


```

```{r functions, include=FALSE}
# cache fit model as rds file (like in brms)
fit <- function(model, file, data, ...) {
  if(file.exists(file)) {
    file %>% 
      readRDS() %>% 
      return()
  } else {
    f <- model %>% 
      rstan::sampling(data = data, ...) 
    
    f %>% saveRDS(file)
    
    return(f)
  }
}

# add energy diagnostic
get_diagnostics <- function(fit) {
  fit %>%
    get_sampler_params(inc_warmup = FALSE) %>%
    map_dfr(as_tibble, .id = 'chain') %>%
    mutate(chain = chain %>% as.integer) %>% 
    group_by(chain) %>% 
    mutate(iter = row_number()) %>%
    ungroup() %>%
    gather(parameter, value, energy__) %>% 
    select(parameter, value, chain, iter) %>% 
    return()
}

# get draws together with energy diagnostic
get_draws <- function(fit, pars) {
  
  if('energy__' %in% pars) {
    energy <- get_diagnostics(fit)
  } else {
    energy <- tibble(
      parameter = character(),
      chain = integer(),
      iter = integer(),
      value = double()
    )
  }
  
  pars <- pars[!str_detect(pars, "energy")]
  
  fit %>% 
    rstan::extract(
      pars,
      permuted=FALSE # to align with energy diagnostic
    ) %>% 
    apply(
      3, 
      as_tibble, 
      .name_repair = ~str_extract(., "\\d+")
    ) %>% 
    map(~mutate(., iter = row_number())) %>%
    bind_rows(.id = 'parameter') %>% 
    gather(chain, value, -parameter, -iter) %>% 
    mutate(chain = chain %>% as.integer()) %>% 
    bind_rows(energy)
}

# add names and indexes of parameters
name_parameters <- function(draws) {
  draws %>% 
    mutate(
      idxs = str_match(parameter, "(\\d+(,\\d+)?)")[,2], 
      name = str_extract(parameter, "^[^\\[]+")
    ) %>% 
    separate(idxs, c('id', 'idx')) 
}
```


## Data Generating Process

Let's recap on the story from last time. We have a 2-year old company that has grown linearly over that time to gain a total of 1000 customers.

```{r customers}
set.seed(65130) # https://www.random.org/integers/?num=2&min=1&max=100000&col=5&base=10&format=html&rnd=new

customers <- tibble(id = 1:1000) %>% 
  mutate(
    end = 2 * 365,
    start = runif(n(), 0, end - 1),
    T = end - start
  )
```

Within a customer's lifetime $\tau$, they will purchase with Poisson-rate $\lambda$. We can simulate the time $t$ till last observed purchase and number of purchases $k$ with `sample_conditional`.

```{r likelihood}
sample_conditional <- function(T, tau, lambda) {
  
  # start with 0 purchases
  t <- 0
  k <- 0
  
  # simulate time till next purchase
  wait <- rexp(1, lambda)
  
  # keep purchasing till end of life/observation time
  while(t + wait <= pmin(T, tau)) {
    t <- t + wait
    k <- k + 1
    wait <- rexp(1, lambda)
  }
  
  # return tabular data
  tibble(
    t = t,
    k = k
  )
}

s <- sample_conditional(300, 200, 1) 

```

```{r sample_conditional, echo=FALSE}
s %>% 
  kable(caption = 'Example output from sample_conditional') %>% 
  kable_styling(bootstrap_options = c("responsive"))
```

In the above example, even though the observation time is $T = 300$, the time $t$ till last purchase will always be below the lifetime $\tau = 200$. With a purchase rate of 1 per unit time, we expect around $k = 200$ purchases.

## Model

We'll use the same likelihood as before, which says that the probability of customer $i$'s data given their parameters is 

$$
\begin{align}
  \mathbb P(k, t, T \mid \mu_i, \lambda_i)
  &=
  \frac{\lambda_i^k}{\lambda_i + \mu_i}
  \left( \mu_i e^{-t(\lambda_i + \mu_i)} + \lambda_i e^{-T(\lambda_i + \mu_i)} \right)
  \\
  &\propto
  p \dpoisson(k \mid t\lambda_i)S(t \mid \mu_i) 
  \\
  &\hphantom{\propto}
  + (1 - p) \dpoisson(k \mid t\lambda_i)\dpoisson(0 \mid (T-t)\lambda_i)S(T \mid \mu_i)
  ,
  \\
  p
  &:=
  \frac{\mu_i}{\lambda_i + \mu_i}
  ,
\end{align}
$$

where $S$ is the exponential survival function.

To turn this into a Bayesian model, we'll need priors for the parameters. The last time, we put simple gamma priors on the parameters $\mu_i$ and $\lambda_i$. For example, we could choose $\lambda_i \sim \dgamma(2, 28)$ if we were to use the simple model from last time (similarly for $\mu_i$). This time we're going hierarchical. There are various ways to make this hierarchical. Let's look at two of them.

[One method](./models/rf.stan) arises directly from the difficulty of specifying the gamma-prior parameters. It involves just turning those parameters into random variables to be simultaneously estimated along with $\mu_i$ and $\lambda_i$. For example, we say $\lambda_i \sim \dgamma(\alpha, \beta)$, where $\alpha_i \sim \dgamma(\alpha_\alpha, \beta_\alpha)$ and $\beta_i \sim \dgamma(\alpha_\beta, \beta_\beta)$ (and similarly for $\mu_i$). We eventually want to incorporate covariates, which is difficult with this parameterisation, so let's move onto a different idea.

Another solution is to use log-normal priors. This means setting 

$$
\begin{align}
  \lambda_i 
  &:= 
  \exp(\alpha_i)
  \\
  \alpha_i
  &\sim
  \dnormal(\beta, \sigma)
  \\
  \beta
  &\sim
  \dnormal(m, s_m)
  \\
  \sigma
  &\sim
  \dnormal_+(0, s)
,
\end{align}
$$ 

where $m$, $s_m$, and $s$ are constants specified by the user (and similarly for $\mu_i$). This implies

* there is an overall mean value $e^\beta$, 
* the customer-level effects $\alpha_i$ are deviations from the overall mean, and
* the extent of these deviations is controlled by the magnitude of $\sigma$.

With $\sigma \approx 0$, there can be very little deviation from the mean, so most customers would be the same. On the other hand, large values of $\sigma$ allow for customers to be (almost) completely unrelated to each other. This means that $\sigma$ is helping us to regularise the model.

The above parameterisation is called "centred", which basically means the prior for $\alpha_i$ is expressed in terms of other parameters ($\beta$, $\sigma$). This can be rewritten as a "non-centred" parameterisation as

$$
\begin{align}
  \lambda_i 
  &:= 
  \exp(\beta + \sigma \alpha_i)
  \\
  \alpha_i
  &\sim
  \dnormal(0, 1)
  \\
  \beta
  &\sim
  \dnormal(m, s_m)
  \\
  \sigma
  &\sim
  \dnormal_+(0, s)
.
\end{align}
$$ 

Notice the priors now contain no references to any other parameters. This is equivalent to the centred parameterisation because $\beta + \sigma \alpha_i \sim \dnormal(\beta, \sigma)$. The non-centred parameterisation is interesting because it is known to increase the sampling efficiency of HMC-based samplers (such as Stan's) in some cases.

## Centred Stan implementation

Here is a [centred stan implementation](./models/rf.stan) of our log-normal hierarchical model. 

```{r rf_centred, results='hide'}
centred <- here::here('models/rf_centred.stan') %>% 
  stan_model()
```

Note that we have introduced the `prior_only` flag. When we specify that we want `prior_only`, then stan will not consider the likelihood and will instead just draw from the priors. This allows us to make prior-predictive simulations. We'll generate a dataset using the prior-predictive distribution, then fit our model to that dataset. The least we can expect from a model is that it fits well to data drawn from its prior distribution. 

### Simulate the dataset

To simulate datasets we'll use hyperpriors that roughly correspond to the priors from the previous post. In particular, the expected lifetime is around 31 days, and the expected purchase rate around once per fortnight.

```{r data_prior}
data_hyperpriors <- list(
  log_life_mean_mu = log(31),
  log_life_mean_sigma = 0.7,
  log_life_scale_sigma = 0.8,

  log_lambda_mean_mu = log(1 / 14),
  log_lambda_mean_sigma = 0.3,
  log_lambda_scale_sigma = 0.5
)

data_prior <- customers %>% 
  mutate(t = 0, k = 0) %>% 
  tidybayes::compose_data(data_hyperpriors, prior_only = 1)

data_prior %>% str()
```

Let's simulate 8 possible datasets from our priors. Notice how the centres and spreads of the datasets can vary.

```{r prior_predictive, results='hide'}
centred_prior <- centred %>% 
  fit( # a wrapper around rstan::sampling to allow caching
    file = here::here('models/rf_centred_prior.rds'), # cache
    data = data_prior,
    pars = c('customer'), # ignore this parameter
    include = FALSE,
    chains = 8,
    cores = 4,
    warmup = 1000, # not sure why this needs to be so high
    iter = 1001, # one more than warmup because we just want one dataset per chain
    seed = 3901 # for reproducibility
  ) 

centred_prior_draws <- centred_prior %>% 
  get_draws( # rstan::extract but also with energy
    pars = c(
    'lp__', 'energy__',
    'theta',
    'log_centres',
    'scales'
    )
  ) %>% 
  name_parameters() # add customer id, and idx = 1 (mu) or 2 (lambda)

```

```{r prior_predictive_distributions, echo=FALSE, fig.cap="Some prior-predictive draws.", fig.height=8, fig.width=8}
centred_prior_draws %>% 
  filter(!is.na(idx)) %>% 
  mutate(
    value = if_else(idx == 1, 1 / value, value),
    name = if_else(idx == 1, '1 / mu', 'lambda')
  ) %>%
  filter(idx == 1 & value < 500 | value < 0.4) %>%
  ggplot() +
  aes(value) +
  geom_histogram() +
  # scale_x_continuous(limits = c(NA, 0.3)) +
  facet_grid(chain ~ name, scales = 'free_x') +
  labs(
    x = 'Value',
    y = 'Count',
    title = 'Histograms of customer-level parameters',
    subtitle = 'for 8 prior-predictive draws',
    caption = "the x-axes have been truncated for clarity"
  )
```

Here are the exact hyperparameters used.

```{r hyper}
hyper <- centred_prior_draws %>% 
  filter(str_detect(parameter, "^log_|scales")) %>% 
  select(chain, parameter, value) %>% 
  spread(parameter, value) 
```

```{r hyper_table, echo=FALSE}
hyper %>%
  kable(caption = "Hyperparameters of the prior-predictive draws.") %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "responsive"))
```

We'll add the prior predictive parameter draws from chain 1 to our customers dataset.

```{r df}
set.seed(33194)

df <- centred_prior_draws %>% 
  filter(chain == 1) %>% 
  filter(name == 'theta') %>% 
  transmute(
    id = id %>% as.integer(), 
    parameter = if_else(idx == '1', 'mu', 'lambda'),
    value
  ) %>% 
  spread(parameter, value) %>% 
  mutate(tau = rexp(n(), mu)) %>% 
  inner_join(customers, by = 'id') %>% 
  group_by(id) %>% 
  group_map(~sample_conditional(.$T, .$tau, .$lambda) %>% bind_cols(.x))

data_df <- data_hyperpriors %>% 
  tidybayes::compose_data(df, prior_only = 0)

```


```{r df_table, echo=FALSE}
df %>% 
  head() %>% 
  kable(caption = 'Sample of customers and their properties') %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "responsive"))
  
```


### Fit the model to simulations

Now we can fit the model to the prior-predictive data `df`.

```{r rf_centred_fit, message=TRUE, warning=TRUE}
centred_fit <- centred %>% 
  fit( # like rstan::sampling but with file-caching as in brms
    file = here::here('models/rf_centred_fit.rds'), # cache
    data = data_df,
    chains = 4,
    cores = 4,
    warmup = 2000,
    iter = 3000,
    control = list(max_treedepth = 12),
    seed = 24207,
    pars = c('customer'),
    include = FALSE
  ) 

centred_fit %>% 
  check_hmc_diagnostics()

```

The HMC diagnostics pass. However, in some of the runs not shown here, there were pretty severe problems with the E-BFMI diagnostic (~0.01) and I've yet to figure out exactly which kinds of situations cause these energy problems. Let's check out the pairwise posterior densities of energy with the hyperparameters.

```{r rf_centred_draws, include=FALSE}
centred_draws <- centred_fit %>% 
  get_draws(
    c('lp__', 'energy__',
      'log_centres', 'scales',
      'theta'
    )
  ) %>% 
  name_parameters()
```

```{r centred_pairwise, fig.height=12, fig.width=16, echo=FALSE, fig.cap="Pairwise posterior densities of the centred model"}
centred_draws %>% 
  filter(parameter == 'energy__' | parameter == 'lp__' | str_detect(parameter, '^(scales|log_centres)')) %>% 
  inner_join(., ., by = c('chain', 'iter')) %>% 
  filter(parameter.x < parameter.y) %>% 
  ggplot() +
  aes(value.x, value.y) +
  geom_hex(aes(colour = ..ndensity.., fill = ..ndensity..)) +
  facet_wrap(~ str_glue("{parameter.y} vs {parameter.x}"), scales = 'free') +
  labs(
    x = 'Right parameter',
    y = 'Left parameter',
    title = 'Pairwise posterior densities',
    subtitle = 'of energy, lp, and hyperparameters',
    fill = 'Density'
  ) +
  guides(colour = FALSE) +
  NULL

```

The scale parameter for the expected lifetime (`scales[1]`) is correlated with energy, which is associated with the energy problems described above. I'm not sure how much of a problem this poses, so let's check out some more diagnostics.

```{r eff_centred}
neff <- centred_fit %>% 
  neff_ratio() %>% 
  tibble(
    ratio = .,
    parameter = names(.)
  ) %>% 
  filter(ratio < 0.5) %>% 
  arrange(ratio) %>% 
  head(20) 
```

```{r eff_centred_table, echo=FALSE} 
neff %>% 
  kable(caption = 'Parameters with the lowest effective sample size.') %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "responsive"))
  
```

Both the `lp__` and `scales[1]` parameters have low effective sample sizes. The rhat values seem fine though.

```{r centred_rhat}
centred_rhat <- centred_fit %>% 
  rhat() %>% 
  tibble(
    rhat = .,
    parameter = names(.)
  ) %>% 
  summarise(
    min_rhat = min(rhat, na.rm = TRUE),
    max_rhat = max(rhat, na.rm = TRUE)
  ) 
```

```{r centred_rhat_table, echo=FALSE}
centred_rhat %>% 
  kable(caption = 'Most extreme rhat values') %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "responsive"))
```

Let's now compare the 90% posterior intervals with the true values. Ideally close to 90% of the 90% posterior intervals capture their true value.

```{r cis_centred}
centred_cis <- centred_draws %>% 
  group_by(parameter) %>% 
  summarise(
    lo = quantile(value, 0.05),
    point = quantile(value, 0.50),
    hi = quantile(value, 0.95)
  ) %>% 
  filter(!str_detect(parameter, '__')) # exclude diagostic parameters
```

The table below shows we managed to recover three of the hyperparameters. The `scales[2]` parameter was estimated slightly too high.

```{r calibration_hyper_centred}
calibration_hyper <- hyper %>% 
  filter(chain == 1) %>% 
  gather(parameter, value, -chain) %>% 
  inner_join(centred_cis, by = 'parameter') %>% 
  mutate(hit = lo <= value & value <= hi)
```

```{r calibration_hyper_centred_table, echo=FALSE}
calibration_hyper %>% 
  select(-chain) %>% 
  kable(caption = "The true hyperparameters and their 90% posterior intervals.") %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "responsive"))

```

We get fairly close to 90% of the customer-level parameters.

```{r, warning=FALSE}
true_values <- df %>% 
  select(id, mu, lambda) %>% 
  gather(parameter, value, -id) %>% 
  mutate(
    idx = if_else(parameter == 'mu', 1, 2),
    parameter = str_glue("theta[{id},{idx}]")
  )

centred_calibration <- centred_cis %>% 
  inner_join(true_values, by = 'parameter') %>% 
  ungroup() %>% 
  summarise(mean(lo <= value & value <= hi)) %>% 
  pull() %>% 
  percent()

centred_calibration
```

This is slightly higher than the ideal value of 90%.

## Non-centred Stan implementation

Here is a [non-centred stan implementation](./models/rf_noncentred.stan) of our log-normal hierarchical model. The important difference is in the expression for $\theta$ and in the prior for $\text{customer}$.

```{r rf_noncentred, results='hide'}
noncentred <- here::here('models/rf_noncentred.stan') %>% 
  stan_model()
```

Since the non-centred and centred models are equivalent, we can also consider `df` as a draw from the non-centred prior predictive distribution. 

```{r rf_noncentred_fit, message=TRUE, warning=TRUE}
noncentred_fit <- noncentred %>% 
  fit( # like rstan::sampling but with file-caching as in brms
    file = here::here('models/rf_noncentred_fit.rds'), # cache
    data = data_df,
    chains = 4,
    cores = 4,
    warmup = 2000,
    iter = 3000,
    control = list(max_treedepth = 12),
    seed = 1259,
    pars = c('customer'),
    include = FALSE
  ) 

noncentred_fit %>% 
  check_hmc_diagnostics()

```

Again, the HMC diagnostics indicate no problems. Let's check the pairwise densities anyway.

```{r rf_noncentred_draws, include=FALSE}
noncentred_draws <- noncentred_fit %>% 
  get_draws(
    c('lp__', 'energy__',
      'log_centres', 'scales',
      'theta'
    )
  ) %>% 
  name_parameters()
```

```{r noncentred_pairwise, fig.height=12, fig.width=16, echo=FALSE, fig.cap="Pairwise posterior densities of the non-centred model"}
noncentred_draws %>% 
  filter(parameter == 'energy__' | parameter == 'lp__' | str_detect(parameter, '^(scales|log_centres)')) %>% 
  inner_join(., ., by = c('chain', 'iter')) %>% 
  filter(parameter.x < parameter.y) %>% 
  ggplot() +
  aes(value.x, value.y) +
  geom_hex(aes(colour = ..ndensity.., fill = ..ndensity..)) +
  facet_wrap(~ str_glue("{parameter.y} vs {parameter.x}"), scales = 'free') +
  labs(
    x = 'Right parameter',
    y = 'Left parameter',
    title = 'Pairwise posterior densities',
    subtitle = 'of energy, lp, and hyperparameters',
    fill = 'Density'
  ) +
  guides(colour = FALSE) +
  NULL

```

The correlation between `scales[1]` and `energy__` is smaller with the non-centred parameterisation. This is reflected in the higher effective sample size for `scales[1]` below. Unfortunately, the effective sample size for the purchase rate hyperpriors has gone down.

```{r eff_noncentred}
neff <- noncentred_fit %>% 
  neff_ratio() %>% 
  tibble(
    ratio = .,
    parameter = names(.)
  ) %>% 
  filter(ratio < 0.5) %>% 
  arrange(ratio) %>% 
  head(20) 
```

```{r eff_noncentred_table, echo=FALSE} 
neff %>% 
  kable(caption = 'Parameters with the lowest effective sample size.') %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "responsive"))
  
```

Again, the rhat values seem fine.

```{r noncentred_rhat}
noncentred_rhat <- noncentred_fit %>% 
  rhat() %>% 
  tibble(
    rhat = .,
    parameter = names(.)
  ) %>% 
  summarise(
    min_rhat = min(rhat, na.rm = TRUE),
    max_rhat = max(rhat, na.rm = TRUE)
  ) 
```

```{r noncentred_rhat_table, echo=FALSE}
noncentred_rhat %>% 
  kable(caption = 'Most extreme rhat values') %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "responsive"))
```


Let's check how many of the 90% posterior intervals contain the true value. 

```{r cis_noncentred}
noncentred_cis <- noncentred_draws %>% 
  group_by(parameter) %>% 
  summarise(
    lo = quantile(value, 0.05),
    point = quantile(value, 0.50),
    hi = quantile(value, 0.95)
  ) %>% 
  filter(!str_detect(parameter, '__')) 
```

The hyperparameter estimates are much the same as with the centred parameterisation.

```{r calibration_hyper_noncentred}
noncentred_calibration_hyper <- hyper %>% 
  filter(chain == 1) %>% 
  gather(parameter, value, -chain) %>% 
  inner_join(noncentred_cis, by = 'parameter') %>% 
  mutate(hit = lo <= value & value <= hi)
```

```{r calibration_hyper_noncentred_table, echo=FALSE}
noncentred_calibration_hyper %>% 
  select(-chain) %>% 
  kable(caption = "The true hyperparameters and their 90% posterior intervals.") %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "responsive"))

```

About 91% of customer-level posterior intervals contain the true value.

```{r, warning=FALSE}
noncentred_calibration <- noncentred_cis %>% 
  inner_join(true_values, by = 'parameter') %>% 
  summarise(mean(lo <= value & value <= hi)) %>% 
  pull() %>% 
  percent()

noncentred_calibration
```

## Discussion

Both centred and non-centred models performed reasonably well on the dataset considered. The non-centred model showed slightly less correlation between `scales` and `energy__`, suggesting it might be the better one to tackle the low E-BFMI problems. Since we only checked the fit on one prior-predictive draw, it would be a good idea to check out the fit to more draws. Some casual attempts of mine (not shown here) suggest there are situations that cause severe E-BFMI problems. Identifying these situations would be an interesting next step. It would also be great to see how it performs on some of the benchmarked datasets mentioned in the [BTYDPlus](https://github.com/mplatzer/BTYDplus) package.
