Import objfull into PSIS notebook, and look at lp__ from tmbstan #41

src/naomi-simple_psis/orderly.yml

@@ -11,6 +11,7 @@ artefacts:
 
 resources:
   - psis.Rmd
+  - naomi_simple.cpp
 
 packages:
   - dplyr
@@ -31,14 +32,11 @@ depends:
       id: latest(parameter:tmb == TRUE && parameter:random_only == TRUE)
       use:
         depends/tmb.rds: out.rds
+        depends/objfull.rds: objfull.rds
   - naomi-simple_fit:
       id: latest(parameter:aghq == TRUE && parameter:k == 3 && parameter:s == 8)
       use:
         depends/aghq.rds: out.rds
-  - naomi-simple_fit:
-      id: latest(parameter:adam == TRUE)
-      use:
-        depends/adam.rds: out.rds
   - naomi-simple_fit:
       id: latest(parameter:tmbstan == TRUE && parameter:niter > 50000)
       use:

src/naomi-simple_psis/psis.Rmd

@@ -24,9 +24,6 @@ Import these inference results as follows:
 ```{r}
 tmb <- readRDS("depends/tmb.rds")
 aghq <- readRDS("depends/aghq.rds")
-adam <- readRDS("depends/adam.rds")
-adam_time <- adam$time
-adam <- adam$adam
 tmbstan <- readRDS("depends/tmbstan.rds")
 
 depends <- yaml::read_yaml("orderly.yml")$depends
@@ -40,50 +37,39 @@ stopifnot(names(tmb$fit$sample) == names(adam$adam$sample))
 stopifnot(names(tmb$fit$sample) == names(tmbstan$mcmc$sample))
 ```
 
-## Run details {.tabset}
-
-For more information about the conditions under which these results were generated, see:
-
-### `TMB`
+# Pareto-smoothed importance sampling
 
-```{r}
-dependency_details <- function(i) {
-  report_name <- names(depends[[i]])
-  print(paste0("Inference results obtained from ", report_name, " with the query ", depends[[i]][[report_name]]$id))
-  report_id <- orderly::orderly_search(query = depends[[i]][[report_name]]$id, report_name)
-  print(paste0("Obtained report had ID ", report_id, " and was run with the following parameters:"))
-  print(orderly::orderly_info(report_id, report_name)$parameters)
-}
-
-dependency_details(1)
-```
+Suppose we have two sets of samples from the posterior.
+For each sample we are going to want to evaluate the log-likelihood, so that we can calculate the log-likelihood ratios.
 
-### `aghq`
+Although we can extract a `TMB` objective function from `tmb` directly, it would correspond to the Laplace approximation.
+As such we use a new objective function using `naomi_simple.cpp` directly:
 
 ```{r}
-dependency_details(2)
+objfull <- readRDS("depends/objfull.rds")
+objfull$fn(objfull$par)
 ```
 
-### `adam`
+The samples from `tmbstan` have the following values of `lp__`.
+These should be equal to what would be obtained by putting the sample into `objfull$fn` (check this):
 
 ```{r}
-dependency_details(3)
+tmbstan_samples <- rstan::extract(tmbstan$mcmc$stanfit)
+
+data.frame(
+  y = tmbstan_samples$lp__,
+  x = 1:length(tmbstan_samples$lp__)
+) %>%
+  ggplot(aes(x = x, y = y)) +
+    geom_point(alpha = 0.4) +
+    labs(y = "Log-likelihood", x = "Draw number") +
+    theme_minimal() 
 ```
 
-### `tmbstan`
+Need to write a function which reformats `TMB` or `tmbstan` samples to be evaluated using `objfull$fn`:
 
 ```{r}
-tmbstan_details <- dependency_details(4)
-tmbstan_details
-```
-
-# Pareto-smoothed importance sampling
+objfull$par
 
-Suppose we have two sets of samples from the posterior.
-For each sample we are going to want to evaluate the log-likelihood, so that we can calculate the log-likelihood ratios.
-We can extract the `TMB` objective function for the log-likelihood as follows:
-
-```{r}
-tmb$fit$obj$fn()
-#' To write!
+#' TODO
 ```