## ----setup_lies,warning=FALSE-------------------------------------------------
library(EcoEnsemble)
set.seed(1234)

## ----eval = FALSE, message=FALSE, warning=FALSE, results='hide', silent=TRUE----
# cor_pri_st <- rstan::stan_model(model_code = " functions{
#   real priors_cor_hierarchical_beta(matrix ind_st_cor, int N, matrix M){
#     real log_prior = 0;
#     for (i in 1:(N-1)){
#       for (j in (i+1):N){
#         log_prior += beta_lpdf(0.5*(ind_st_cor[i, j] + 1)| M[i, j], M[j, i] );
#       }
#     }
#     return log_prior;
#   }
# 
#   real beta_conj_prior(real alpha, real betas, real r, real s, real k){
#     real rl = 1/(1 + exp(-r));
#     real sl = 1/(1 + exp(-s));
#     real p = (sl * rl)^k;
#     real q = (sl * (1 - rl))^k;
#     real ret = alpha * log(p) + betas * log(q) - k * lbeta(alpha,betas);
#     return ret;
#   }
# 
# }
# 
# data {
#   vector[3] cor_p;
# }
# 
# parameters {
#   matrix <lower=0>[5,5] beta_pars;
#   corr_matrix[5] rho[4];
# }
# 
# model {
#   for (i in 1:4){
#     for (j in (i+1):5){
#       target += beta_conj_prior(beta_pars[i,j], beta_pars[j,i], cor_p[1], cor_p[2], cor_p[3]);
#     }
#   }
# 
#   for (i in 1:4){
#     target += priors_cor_hierarchical_beta(rho[i],5,beta_pars);
#     diagonal(beta_pars) ~ std_normal();
#   }
# 
# }
# 
# generated quantities {
# 
#   matrix [5,5] rhovars;
#   for (i in 1:4){
#     for (j in (i+1):5){
#       rhovars[i,j] = 4 * (beta_pars[i,j] * beta_pars[j,i])/(square(beta_pars[i,j] + beta_pars[j,i]) * (beta_pars[i,j] + beta_pars[j,i] + 1));
#       rhovars[j,i] = (2 * beta_pars[i,j]/(beta_pars[i,j] + beta_pars[j,i])) - 1;
#     }
#   }
# 
#   for (i in 1:5){
#     rhovars[i,i] = 4;
#   }
# 
# }
# ")
# library(ggplot2)
# rhoplots <- list()
# kvals <- c(0.05, 5)
# parvals <- do.call(expand.grid, c(rep(list(c(0.25, 3)), 2), list(kvals)))
# #Sampling and gathering outputs for plotting
# for (i in 1:nrow(parvals)){
#   fit_cor <- rstan::sampling(cor_pri_st, data = list(cor_p=as.numeric(parvals[i,])), iter = 2000, chains=4)
#   ex.fit <- rstan::extract(fit_cor)
#   rho_density <- density(ex.fit$rho[,1,1,2], from = -1, to = 1)
#   rhoplot_data <- data.frame(rho_density$x, rho_density$y)
#   names(rhoplot_data) <- c("rho", "Density")
#   rhoplot_data <- cbind(rhoplot_data, rep(parvals[i,1], nrow(rhoplot_data)), rep(parvals[i,2], nrow(rhoplot_data)), rep(parvals[i,3], nrow(rhoplot_data)))
#   names(rhoplot_data)[3:5] <- c("r", "s", "k")
#   rhoplots[[i]] <- rhoplot_data
# }
# #Construct plots
# rhoplot_data <- do.call(rbind, rhoplots)
# rhoplot_data <- cbind(rhoplot_data, rep("Beta conjugate prior", nrow(rhoplot_data)))
# names(rhoplot_data)[6] <- "Prior"
# unif_range <- seq(-1, 1, length.out = nrow(rhoplots[[1]]))
# unif_data <- data.frame(cbind(rep(unif_range, nrow(parvals)), rep(dbeta((unif_range+1)/2,5/2,5/2)/2, nrow(parvals)), rhoplot_data[,3:5], rep("Uniform prior", nrow(rhoplot_data))))
# names(unif_data) <- names(rhoplot_data)
# rhoplot_data <- rbind(rhoplot_data, unif_data)
# rhoplot1 <- ggplot(rhoplot_data[which(rhoplot_data$k == kvals[1]),]) + geom_area(aes(x = rho, y = Density, color = Prior, fill = Prior), alpha = 0.3, position = "identity") + scale_x_continuous(bquote(rho), sec.axis = sec_axis(~., name = "r", breaks = NULL, labels = NULL)) + scale_y_continuous(sec.axis = sec_axis(~., name = "s", breaks = NULL, labels = NULL)) + theme(aspect.ratio = 1) + facet_grid(rows = vars(r), cols = vars(s)) + scale_color_manual(values = c("blue", "green")) + scale_fill_manual(values = c("blue", "green"))
# rhoplot2 <- ggplot(rhoplot_data[which(rhoplot_data$k == kvals[2]),], aes(x = rho, y = Density)) + geom_area(aes(color = Prior, fill = Prior), alpha = 0.3, position = "identity") + scale_x_continuous(bquote(rho), sec.axis = sec_axis(~., name = "r", breaks = NULL, labels = NULL)) + scale_y_continuous(sec.axis = sec_axis(~., name = "s", breaks = NULL, labels = NULL)) + theme(aspect.ratio = 1) + facet_grid(rows = vars(r), cols = vars(s)) + scale_color_manual(values = c("blue", "green")) + scale_fill_manual(values = c("blue", "green"))

## ----echo = FALSE, out.width="700px", out.height="700px"----------------------
knitr::include_graphics("data/rhoplot1.png")

## ----eval = FALSE, fig.dim = c(7,6), echo = FALSE-----------------------------
# rhoplot1

## ----echo = FALSE, out.width="700px", out.height="700px"----------------------
knitr::include_graphics("data/rhoplot2.png")

## ----eval = FALSE, fig.dim = c(7,6), echo = FALSE-----------------------------
# rhoplot2

## ----message=FALSE, warning=FALSE,eval=FALSE, results='hide', silent=TRUE-----
# kvals <- c(0.1, 0.2, 0.4, 0.8, 1.6, 3.2)
# #parvals <- cbind(rep(0.3, 6), rep(-1, 6), kvals)
# parvals <- cbind(rep(-0.3, 6), rep(3, 6), kvals)
# rhovarplots <- list()
# rhomeanplots <- list()
# for (i in 1:6){
#   fit_cor <- rstan::sampling(cor_pri_st, data = list(cor_p=parvals[i,]), iter = 2000, chains=4)
#   ex.fit <- rstan::extract(fit_cor)
#   rhovar_density <- density(ex.fit$rhovars[,1,2])
#   rhomean_density <- density(ex.fit$rhovars[,2,1])
#   rhovar_data <- data.frame(cbind(rhovar_density$x, rhovar_density$y, rep(kvals[i], length(rhovar_density$x))))
#   rhomean_data <- data.frame(cbind(rhomean_density$x, rhomean_density$y, rep(kvals[i], length(rhomean_density$x))))
#   names(rhovar_data) <- c("Variance", "Density", "k")
#   names(rhomean_data) <- c("Expectation", "Density", "k")
#   rhovarplots[[i]] <- rhovar_data
#   rhomeanplots[[i]] <- rhomean_data
# }
# rhovarplot_data <- do.call(rbind, rhovarplots)
# rhomeanplot_data <- do.call(rbind, rhomeanplots)
# rhovarplot <- ggplot(rhovarplot_data) + geom_area(aes(x = Variance, y = Density), color = "blue", fill = "blue", alpha = 0.3, position = "identity") + geom_vline(xintercept = 0.2, color = "green", linetype = "dashed", linewidth = 0.8) + facet_wrap(vars(k), nrow = 2, ncol = 3) + scale_x_continuous(sec.axis = sec_axis(~., name = "k", breaks = NULL, labels = NULL))
# rhomeanplot <- ggplot(rhomeanplot_data) + geom_area(aes(x = Expectation, y = Density), color = "blue", fill = "blue", alpha = 0.3, position = "identity") + geom_vline(xintercept = 0, color = "green", linetype = "dashed", linewidth = 0.8) + facet_wrap(vars(k), nrow = 2, ncol = 3) + scale_x_continuous(sec.axis = sec_axis(~., name = "k", breaks = NULL, labels = NULL))

## ----echo = FALSE, out.width="600px", out.height = "500px"--------------------
knitr::include_graphics("data/rhovarplot.png")

## ----eval = FALSE, echo = FALSE-----------------------------------------------
# rhovarplot

## ----echo = FALSE, out.width = "600px", out.height = "500px"------------------
knitr::include_graphics("data/rhomeanplot.png")

## ----eval=FALSE, fig.dim = c(8,5), echo = FALSE-------------------------------
# rhomeanplot

## ----eval = FALSE, message = FALSE, warning = FALSE, silent = TRUE, results = 'hide'----
# fit_cor <- rstan::sampling(cor_pri_st, data = list(cor_p=c(0.25, 3, 4), iter = 2000, chains=4))
# ex.fit <- rstan::extract(fit_cor)
# rhoplot_data  <- data.frame(ex.fit$rho[,1,1,2])
# names(rhoplot_data) <- "rho"
# rhoplot <- ggplot(rhoplot_data) + geom_histogram(aes(x = rho), color = "blue", fill = "blue", alpha = 0.3, binwidth = 0.05) + scale_x_continuous(bquote(rho))

## ----echo = FALSE, out.width = "400px", out.height = "400px"------------------
knitr::include_graphics("data/rhoplot.png")

## ----eval = FALSE, echo = FALSE, warning = FALSE, fig.dim = c(4,3)------------
# rhoplot

## ----eval=FALSE,message=FALSE, warning = FALSE, results = 'hide',silent=TRUE----
# priors <- EnsemblePrior(4,
#                         ind_st_params =IndSTPrior("hierarchical_beta_conjugate",list(-3,1,8,4),
#                                                   list(0.25,3,4),AR_params=c(2,2)),
#                         ind_lt_params = IndLTPrior("lkj",list(1,1),1),
#                         sha_st_params = ShaSTPrior("lkj",list(1,10),1,AR_params=c(2,2)),
#                         sha_lt_params = 5
# )
# prior_density <- prior_ensemble_model(priors, M = 4)
# samples <- sample_prior(observations = list(SSB_obs, Sigma_obs),
#              simulators = list(list(SSB_ewe, Sigma_ewe,"EwE"),
#                    list(SSB_fs ,  Sigma_fs,"FishSums"),
#                    list(SSB_lm ,  Sigma_lm,"LeMans"),
#                    list(SSB_miz, Sigma_miz,"mizer")),
#              priors=priors,
#              sam_priors = prior_density)

## ----eval = FALSE, fig.dim = c(7, 6)------------------------------------------
# plot(samples)

## ----echo = FALSE, out.height="500px", out.width="800px"----------------------
knitr::include_graphics("data/p_priors_beta_conj.png")