Given a deviance function \(D(\eta) = -2 \log(p(y|\eta))\), and an estimate
\(\eta* = (\sum \eta_i) / n\) of the posterior mean
\(E(\eta|y)\), where \(y = (y_1, ..., y_n)\) denote the data, \(\eta\) is the unknown
parameter vector of the model, \(\eta_1, ..., \eta_N\) are MCMC samples from the posterior
distribution of \(\eta\) given \(y\) and \(p(y|\eta)\) is the likelihood function,
the Watanabe-Akaike Information Criterion (WAIC) is defined as
$$WAIC = LPPD - p_W$$
where
$$LPPD = \sum_{i=1}^n \log (N^{-1} \sum_{s=1}^N p(y_i|\eta_s) )$$
and (form 1 of)
$$p_W = 2 \sum_{i=1}^n [ \log (N^{-1} \sum_{s=1}^N p(y_i|\eta_s) ) - N^{-1} \sum_{s=1}^N \log \:p(y_i|\eta_s) ].$$
An alternative form (form 2) for \(p_W\) is given by
$$p_W = \sum_{i=1}^n \hat{var} \log p(y_i|\eta)$$
where for \(i = 1, ..., n\), \(\hat{var} \log p(y_i|\eta)\) denotes the estimated variance
of \(\log p(y_i|\eta)\) based on the realizations \(\eta_1, ..., \eta_N\).
Note that waic.angmcmc calls waic for computation. If the likelihood contribution of each data
point for each MCMC iteration is available in object (can be returned by setting return_llik_contri = TRUE)
during fit_angmix call), waic.array is used; otherwise waic.function is
called. Computation is much faster if the likelihood contributions are available - however, they are very
memory intensive, and by default not returned in fit_angmix.