s_plot_results.Rmd

---
title: "Plot electrochemistry"
output: html_notebook
---

```{r}
library(tidyverse)
```

```{r}
get_model_simulations <- function(is_autocorrelated) {
  if(is_autocorrelated)
    filename <- "mcmc_results_autocorrelated.csv"
  else
    filename <- "mcmc_results_independent.csv"
  vals <- read.csv(filename, header = FALSE)$V1
  vals
}

df <- read.csv("data.csv", header = FALSE)
colnames(df) <- c("time", "current")

vals_a <- get_model_simulations(TRUE)
vals_i <- get_model_simulations(FALSE)
df <- df %>% 
  mutate(
    autocorrelated=vals_a,
    independent=vals_i
  )

df %>% 
  filter(time < 20) %>% 
  pivot_longer(-c(time, current)) %>% 
  mutate(name=as.factor(name)) %>% 
  mutate(name=fct_relevel(name, "independent", "autocorrelated")) %>% 
  ggplot(aes(x=time, y=value)) +
  geom_line(aes(y=current)) +
  geom_line(colour="blue") +
  theme_classic() +
  facet_wrap(~name)
```

Examine posteriors
```{r}
get_posteriors <- function(is_autocorrelated) {
  if(is_autocorrelated) {
    filename <- "mcmc_draws_autocorrelated.csv"
    type="autocorrelated"
    params <- c('k0', 'E0', 'a', 'Ru', 'Cdl', 'freq', 'rho', 'sigma')
  } else {
    filename <- "mcmc_draws_independent.csv"
    params <- c('k0', 'E0', 'a', 'Ru', 'Cdl', 'freq', 'sigma')
    type="independent"
  }
  df <- read.csv(filename, header = FALSE)
  colnames(df) <- params
  df %>% 
    mutate(type=type)
}

df_a <- get_posteriors(TRUE)
df_i <- get_posteriors(FALSE)
df_both <- df_i %>% 
  bind_rows(df_a) %>% 
  select(-rho) %>% 
  mutate(iteration=seq_along(type)) %>% 
  pivot_longer(-c(iteration, type))
df_summary <- df_both %>% 
  group_by(type, name) %>% 
  summarise(
    lower=quantile(value, 0.025),
    upper=quantile(value, 0.975),
    middle=quantile(value, 0.5)
  ) %>% 
  pivot_wider(id_cols = name,
              values_from = c(lower, upper, middle),
              names_from = type)

df_summary %>% 
  ggplot(aes(x=middle_independent, y=middle_autocorrelated,
             colour=name)) +
  geom_abline() +
  geom_pointrange(aes(ymin=lower_autocorrelated,
                      ymax=upper_autocorrelated)) +
  geom_errorbarh(aes(xmin=lower_independent,
                     xmax=upper_independent),
                 height=0) +
  scale_x_log10() +
  scale_y_log10() +
  scale_color_brewer("Parameter", palette = "Dark2") +
  xlab("Independent estimate") +
  ylab("Autocorrelated estimate")
```
# Results from a bigger variety of models

```{r}
get_model_simulations <- function(error_model) {
  filename <- paste0("mcmc_results_", error_model, ".csv")
  vals <- read.csv(filename, header = FALSE)$V1
  vals
}

vals_arma11 <- get_model_simulations("arma11")
vals_ar1 <- get_model_simulations("ar1")
vals_i <- get_model_simulations("independent")
vals_s <- get_model_simulations("student_t")

df <- read.csv("data.csv", header = FALSE)
colnames(df) <- c("time", "current")

df <- df %>% 
  mutate(
    independent=vals_i,
    "AR(1)"=vals_ar1,
    "ARMA(1,1)"=vals_arma11,
    "Student-t"=vals_s
  )

g1 <- df %>% 
  filter(time < 25) %>% 
  pivot_longer(-c(time, current)) %>% 
  mutate(name=as.factor(name)) %>% 
  mutate(name=fct_relevel(name, "independent", "AR(1)", "ARMA(1,1)", "Student-t")) %>% 
  ggplot(aes(x=time, y=value)) +
  geom_line(aes(y=current)) +
  geom_line(colour="blue") +
  theme_classic() +
  facet_wrap(~name) +
  xlab("Time") +
  ylab("Current")
g1
```


```{r}
get_posteriors <- function(error_model, params) {
  filename <- paste0("mcmc_draws_", error_model, ".csv")
  df <- read.csv(filename, header = FALSE)
  colnames(df) <- params
  df %>% 
    mutate(type=error_model)
}

params <- c('k0', 'E0', 'a', 'Ru', 'Cdl', 'freq', 'sigma')
df_i <- get_posteriors("independent", params)
params <- c('k0', 'E0', 'a', 'Ru', 'Cdl', 'freq', 'rho', 'sigma')
df_ar1 <- get_posteriors("ar1", params)
params <- c('k0', 'E0', 'a', 'Ru', 'Cdl', 'freq', 'rho', 'phi', 'sigma')
df_arma11 <- get_posteriors("arma11", params)
params <- c('k0', 'E0', 'a', 'Ru', 'Cdl', 'freq', 'nu', 'sigma')
df_s <- get_posteriors("student_t", params)

df_all <- df_i %>% 
  bind_rows(
    df_ar1,
    df_arma11,
    df_s
  ) %>% 
  group_by(type) %>% 
  mutate(iteration=seq_along(type)) %>% 
  ungroup() %>% 
  pivot_longer(-c(iteration, type))

df_summary <- df_all %>% 
  group_by(type, name) %>% 
  summarise(
    lower=quantile(value, 0.01, na.rm=TRUE),
    middle=quantile(value, 0.5, na.rm=TRUE),
    upper=quantile(value, 0.99, na.rm=TRUE)
  )

df_ind_only <- df_summary %>% 
  filter(type=="independent") %>% 
  pivot_wider(id_cols = name,
              values_from = c(lower, middle, upper),
              names_from = type)
df_wo_ind <- df_summary %>% 
  filter(type!="independent")
df_combined <- df_wo_ind %>% 
  left_join(df_ind_only)

g2 <- df_combined %>% 
  filter(!name %in% c("nu", "phi", "rho")) %>% 
  mutate(type=case_when(
    type=="ar1"~"AR(1)",
    type=="arma11"~"ARMA(1, 1)",
    type=="student_t"~"Student-t"
  )) %>% 
  ggplot(aes(x=middle_independent, y=middle,
             colour=name)) +
  geom_abline() +
  geom_pointrange(aes(ymin=lower,
                      ymax=upper), size=0.1) +
  geom_errorbarh(aes(xmin=lower_independent,
                     xmax=upper_independent),
                 height=0) +
  scale_x_log10() +
  scale_y_log10() +
  scale_color_brewer("Parameter", palette = "Dark2") +
  xlab("Independent error estimate") +
  ylab("Alternative error model estimate") +
  facet_grid(cols=vars(type))
g2
```


```{r}
library(cowplot)
gfinal <- plot_grid(g1, g2, nrow = 2, labels = c("A.", "B."))
ggsave("error_model_comparison.pdf", gfinal, width = 10, height = 6)
```

Residuals! Hmm. Not sure what to plot since all models will have autocorrelated residuals...
```{r}
df_acf <- df %>% 
  pivot_longer(-c(time, current)) %>% 
  mutate(res=current-value) %>% 
  group_by(name) %>% 
  summarise(
    ac=acf(res)$acf
  )
```