Odd behaviour for some parameters, e.g.:

Likely doesn't make sense that TMB would have larger standard deviation here as aghq should be taking into account more variation (mixture of Gaussians rather than a single Gaussian).

First check that TMB is actually Gaussian. Think this is done. Obtained the joint Hessian using this. Calling to naomi:::sdreport_joint_precision. Looks like no bias correction used etc.

After optimisation, compution of mean vector and Hessian for each grid point is embarrisingly parallel

• Scope out the possibility for more control over the evaluation of MakeADFun, especially the matrix algebra. This can be started by first understanding the internals of MakeADFun. Could it be replaced with relying on and and the Laplace approximation computed outside of TMB, with parts replaced to speed it up as required
  • Put MakeADFun into development script
  • Click through MakeADFun and write comments as to what is doing which part. Try to find the part which is doing the matrix algebra.

• Test how much longer AGHQ hyperparameter quadrature takes as a function of number of grid points
• Test how much longer a given Laplace marginal takes as a function of number of grid points
• Test whether accuracy of latent field Laplace approximations increases
• Test whether accuracy of mixture of Gaussians increases

(As opposed to hyperparamters or latent field nodes)

Jeff points to a repo of Kasper's looking for "Simple examples for which the Laplace approximation is inaccurate":

It links to this git repo with an example of a false convergence error and unpacking it a bit: https://github.com/kaskr/laplace_accuracy

The model that Kasper suggests is causing trouble is:

u1 ~ N(mu1, sd1)
u2 ~ N(mu2, sd2)
N ~ Pois(exp(u1) + exp(u2))

Should try this out with aghq and try to develop intuition as to why it might do better than TMB in this case.

• https://bayescomp2023.com/
• Finished by 12th March
• docs_bayescomp-poster

• Epidemiology e.g. compartmental model
• Seth's work
• Existing TMB implementation?
• Look at TMB examples

• Add KS statistic for all indices i in 1:30
• Plot KS of the Gaussian against KS of the Laplace with an x = y line for reference
• Check that tmbstan isn't sampling from h even though it's passed h not f
• Resolve bug causing notebook i = 1 to be different
• Move clean version to script
• Extend away from empirical Bayes to integration with aghq over the hyperparameters
• Make distributions plot larger in notebook
• Write-up maths into paper.Rmd

• Future development of risk category model likely going to want to be more flexible
• Would be a good case study for inference methods work

• Read Bayesian workflow for disease transmission modeling in Stan (code here) and copy their Stan model to TMB
• Is a compartmental model an ELGM? Probably not right?
• Tim uses https://en.wikipedia.org/wiki/Euler_method as standard in HIV applications since it's fast and performs fine since things change really slowly over time. Maybe it's less suitable for other applications

e.g.

• k = 5 and s = 5 would be 3125 points and take ~0.5 hours
• k = 5 and s = 6 would be 3125 points and take ~1.5 hours

How do these do in the KS?

• Which parameters are poorly identified?
• Are the posteriors from TMB systematically wider than those from aghq or tmbstan: perhaps they may be because less uncertainty in hyperparamters is taken into account

• Make reports to generate these?

This is for both TMB and aghq. What's going on here?

Obtain the joint Hessian using this. Calling to `naomi:::sdreport_joint_precision. Looks like no bias correction used etc.

• Using Simplified INLA (Wood)
• To be submitted as pull request to aghq probably
• Allows testing of replicated INLA

From Alex:

is there anything that stands out to you

I think looking at the mean and SD is a bit too simplistic. With large data and under "standard" conditions, the mean and the mode will be similar (think Bernstein von-Mises, convergence of the posterior to a Gaussian), so I would argue that any approximation method really should be able to pick up the mean well. And, being interested at all the in the second moment sort of implies you believe the posterior is Gaussian anyways, in which case aghq would be exact. So it's definitely encouraging that the two approximations (aghq and MCMC) are close in the first two moments, but it doesn't really tell the whole story.

I think it would be more instructive to measure something to do with the tails. I usually compare the KS distance (maximal difference in empirical CDF) between the samples returned by tmbstan, and those from aghq::sample_marginal, marginally for each parameter. You can't fool the KS– if the tails aren't captured, it will be large. You can use ks.test() in R which also gives a p-value, which irrespective of philosophy and distributional assumptions, does give you a rough measure of what KS values are "large". It's also directly interpretable, like KS = 0.05 means that no matter where you compute a tail probability, you're never more than 5 percentage points away from what you'd get using tmbstan. I do this in the software vignette, https://arxiv.org/pdf/2101.04468.pdf, Table 3 for example.

The KS stats only capture the marginal distributions. For the joint, I haven't yet explored this too much, but one potentially compelling option is Pareto-Smoothed Importance Sampling (PSIS, see for example https://arxiv.org/pdf/1802.02538.pdf). You can use the loo::psis function for this, I have gotten it to work before but I can't seem to find the code (sorry about that). It gives you an interpretable, scalar summary of the agreement between the approximate and true joint distributions, which is a very useful thing especially in high dimensions.

The ESS is just be a function of the chains, so should be able to calculate this externally to tmbstan.

Need to fix whatever is causing this! Pretty clear that it's incorrect.

If going to use tmbstan as gold standard, need to check that things look alright (and eventually use longer chains)

Does EB work fine? Could also try Laplace marginals here.

For example, which latent field elements generate a given output?

• Expose internals of fit_aghq and run through
• Create TMB template and objective function for one named latent field parameter and index
• Attain Laplace marginal for one named latent field parameter with empirical Bayes hypers
• Attain Laplace marginal for one named latent field parameter with grid hypers
• Generalise code to sweep over all indices of one named latent field parameter
  • I decided this doesn't need to be done, and we can just sweep over all indices once all the named latent field parameters are together
• Generate marginals for beta_rho[1] to compare to -- just use those from naomi-simple_compare
• Understand what's going wrong when the laplace_marginal errors due to the negative weighted posterior outweighing the positive weighted posterior in some location
  • Is there any way to sample from the mixture of Gaussians with a sparse grid and possibly negative weights?
• Start the optimisation within for(z in ...) loop for each evaluation of obj$fn as as close to the eventual optima as possible -- likely at the mode of the Gaussian approximation -- by using the random.start option
  • Benchmark exactly how much cheaper that has made things
• Use the suggested cheaper code to calculate the diagonal of the inverse Hessian
• Generalise code to sweep over all latent field parameters
  • Strategy to do this is to concatenate all the names into x_i and x_minus_i, then reconstruct x, then from x generate all the named parameters based on their position
• Create function fit_name_here which brings this all together
• Fit k = 1 with Laplace marginals
• Fit k = 2 or k = 3 sparse with Laplace marginals
  • I think the sparsity causes a problem for the weights
• A way to sample from the spline interpolated Laplace marginals
• Create sample_adam function to sample from i = 1, ... spline interpolated Laplace marginals
• Compare output from sample_adam to the other inference methods e.g. using KS tests
• #23
• Expose internals of sample_aghq to get a feel again for how the latent field and hyperparameter marginals are being used and run through
  • Blocker: can't fit product grid aghq with k > 1.
    • Possible solution: lower area_level to 3 rather than 4. But this won't fix exponential in number of hypers problem
    • Possible solution: create version of naomi_simple.cpp with fewer hypers. But this would be annoying to do
    • Possible solution: find a way to fix some of the hypers, probably with the map option
• See if it's possible to use DATA_UPDATE macros to avoid recreating TMB template for each value i
  • This looks to be impossible. TMB doesn't like that I if-else branch after updating it:
    
• #22

Howes, Adam T
Does someone know if it's possible to generate data from the prior with tmbstan? Something like doing this to get an objective function without any contribution from the likelihood (but data can't beNA):

Wolock, Tim
I don't know of a built-in way in TMB or Stan. You really just want "y", not "X", to be NA, right? Could you add a use_likelihood object to your data list? You would have something like if (use_likelihood==1) nll -= dnorm(y, mu, sigma, true);. If you "fit" the model, you should just get prior samples of the distributional params, right? If you're generating posterior predictive samples in a SIMULATE block, you could also get prior predictive samples from TMB with obj$simulate(x). (Taking this opportunity to plug numpyro, which makes conditioning optional)

Howes, Adam T
Thanks Tim! Yeah I think that could work about adding a flag. Then whatever the values of y you input it doesn't matter, they count for nothing. I also didn't know about the SIMULATE block, so helpful to know. For some reason with Stan I thought it was easy to generate data from the prior, but I think I was mistaken about this.

Probably better to be integrating over the ECDF difference than finding the maximum.

https://en.wikipedia.org/wiki/Anderson%E2%80%93Darling_test

Could help to determine normalising constant stabilisation.