Update paper analysis

geco-bern · Dec 2, 2024 · 36f990b · 36f990b
1 parent a09b439
commit 36f990b
Show file tree

Hide file tree

Showing 2 changed files with 18 additions and 66 deletions.
diff --git a/analysis/02-bayesian-calibration.R b/analysis/02-bayesian-calibration.R
@@ -60,12 +60,12 @@ plot_prior_posterior_density <- function(
     ggplot(
       aes(x = value, fill = distrib)
     ) +
-    geom_density() +
+    geom_density(adjust = 5) +
     theme_classic() +
     facet_wrap( ~ variable , nrow = 2, scales = "free") +
     theme(legend.position = "bottom",
           axis.title.x = element_text("")) +
-    scale_fill_manual(values = c("#29a274ff", t_col("#777055ff"))) # GECO colors
+    scale_fill_manual(NULL, values = c("#29a274ff", t_col("#777055ff"))) # GECO colors
 
   return(gg)
 }
@@ -80,8 +80,8 @@ settings_calib <- list(
   control = list(
     sampler = "DEzs",
     settings = list(
-      burnin = 3000,
-      iterations = 12000,
+      burnin = 6000,
+      iterations = 18000,
       nrChains = 3,       # number of independent chains
       startValue = 3      # number of internal chains to be sampled
     )),
@@ -114,6 +114,7 @@ toc() # Stop measuring time
 
 # Plot prior and posterior distributions
 gg <- plot_prior_posterior_density(par_calib$mod)
+gg
 
 get_runtime <- function(par_calib) {# function(settings_calib){
   total_time_secs <- sum(unlist(lapply(

diff --git a/analysis/03-uncertainty-estimation.R b/analysis/03-uncertainty-estimation.R
@@ -23,7 +23,7 @@ set.seed(2023)
 
 # Sample parameter values from the posterior distribution
 samples_par <- getSample(par_calib$mod,
-                         thin = 46,              # get 600 samples in total
+                         thin = 60,              # get 600 samples in total
                          whichParameters = 1:4) |>
   as.data.frame() |>
   dplyr::mutate(mcmc_id = 1:n()) |>
@@ -51,7 +51,7 @@ run_pmodel <- function(sample_par){
   # return modelled GPP and prediction for a new GPP observation
   gpp <- out$data[[1]][, "gpp"]
   data.frame(gpp = gpp, 
-             gpp_pred = gpp + rnorm(n = length(gpp), mean = 0, 
+             gpp_pred = rnorm(n = length(gpp), mean = gpp,
                                     sd = sample_par$err_gpp),
              date = out$data[[1]][, "date"])
 }
@@ -88,15 +88,13 @@ plot_gpp_error <- ggplot(data = data_to_plot) +
         ymax = gpp_pred_q95,
         x = date,
         fill = "Model uncertainty"
-    ),
-    alpha = 0.5) +
+    )) +
   geom_ribbon(
     aes(ymin = gpp_q05, 
         ymax = gpp_q95,
         x = date,
         fill = "Parameter uncertainty"
-    ),
-    alpha = 0.5) +
+    )) +
   # Include observations in the plot
   geom_point(
     aes(x = date,
@@ -108,8 +106,7 @@ plot_gpp_error <- ggplot(data = data_to_plot) +
     aes(x = date,
         y = gpp_q50,
         color = "Predictions"
-    ),
-    alpha = 0.6
+    )
   ) +
   theme_classic() +
   theme(panel.grid.major.y = element_line(),
@@ -120,64 +117,18 @@ plot_gpp_error <- ggplot(data = data_to_plot) +
   )
 
 plot_gpp_error <- plot_gpp_error +
-  scale_color_manual(name = "",
+  scale_color_manual(NULL,
                      breaks = c("Observations",
                                 "Predictions"),
                      values = c(t_col("black", 10),
-                                t_col("grey40", 10)))# +
-#  scale_fill_manual(name = "",
-#                    breaks = c("Model uncertainty",
-#                               "Parameter uncertainty"),
-#                    values = c(t_col("#E69F00", 60),
-#                               t_col("#009E73", 40)))
+                                t_col("#7570b3", 10))) +
+  scale_fill_manual(NULL,
+                    breaks = c("Model uncertainty",
+                               "Parameter uncertainty"),
+                    values = c(t_col("#d95f02", 50),
+                               t_col("#1b9e77", 50)))
 plot_gpp_error
+
 settings_string <- get_settings_str(par_calib)
 ggsave(paste0("./analysis/paper_results_files/",settings_string,"_gpp_predictions_observations.pdf"), plot = plot_gpp_error, width = 6, height = 5)
 ggsave(paste0("./analysis/paper_results_files/",settings_string,"_gpp_predictions_observations.png"), plot = plot_gpp_error, width = 6, height = 5)
-
-
-plot_gpp_error2 <- ggplot(data = data_to_plot) +
-  # geom_ribbon(
-  #   aes(ymin = gpp_pred_q05,
-  #       ymax = gpp_pred_q95,
-  #       x = date,
-  #       fill = "Model uncertainty")) +
-  geom_ribbon(
-    aes(ymin = gpp_q05, 
-        ymax = gpp_q95,
-        x = date,
-        fill = "Parameter uncertainty")) +
-  # geom_line(
-  #   aes(x = date,
-  #       y = gpp_q50,
-  #       color = "Predictions")) +
-  # Include observations in the plot
-  # geom_point(
-  #   aes(x = date,
-  #       y = gpp_obs,
-  #       color = "Observations")) +
-  theme_classic() +
-  theme(panel.grid.major.y = element_line(),
-        legend.position = "bottom") +
-  labs(
-    x = 'Date',
-    y = expression(paste("GPP (g C m"^-2, "s"^-1, ")"))
-  )
-
-plot_gpp_error2 <- plot_gpp_error2 +  
-  scale_color_manual(name = "",
-                     breaks = c("Observations",
-                                "Predictions",
-                                "Non-calibrated predictions"),
-                     values = c(t_col("black", 10),
-                                t_col("#E69F00", 10),
-                                t_col("#56B4E9", 10))) +
-  scale_fill_manual(name = "",
-                    breaks = c("Model uncertainty",
-                               "Parameter uncertainty"),
-                    values = c(t_col("#E69F00", 60),
-                               t_col("#009E73", 40)))
-plot_gpp_error2
-ggsave(paste0("./analysis/paper_results_files/",settings_string,"_gpp_predictions_observations_2.pdf"), plot = plot_gpp_error2, width = 6, height = 5)
-ggsave(paste0("./analysis/paper_results_files/",settings_string,"_gpp_predictions_observations_2.png"), plot = plot_gpp_error2, width = 6, height = 5)
-