diff --git a/Fish 558 Homeworks/HWK4/Homework 4 Functions and Data/FISH 558 HWK4.Rmd b/Fish 558 Homeworks/HWK4/Homework 4 Functions and Data/FISH 558 HWK4.Rmd
diff --git a/Fish 558 Homeworks/HWK4/Homework 4 Functions and Data/mozzy_likelihood.R b/Fish 558 Homeworks/HWK4/Homework 4 Functions and Data/mozzy_likelihood.R
@@ -0,0 +1,54 @@
+#' A function to return the likelihood for assignmen
+#'
+#' \code{mosquito.like} returns the negative log likelood of the data
+#' @param pars is a list of parameters for the model
+#' @param dat the stream-distance-count data
+#' @param parnames vector of parameter names
+#' @param stream_pars index of parameters that are stream specific
+#' @param dat data
+#' @param sigma.hyperprior.mu
+
+mozzy_like <- function(pars ,parnames, par_mat,stream_pars, dat , sigma.hyperprior.mu = 1,sigma.hyperprior.sd = 0.1,
+                       alpha.hyperprior.mu = .01, alpha.hyperprior.sd = 0.1)
+{
+
+  # Data likelihood ----
+
+  pred_dat <- dat
+
+  pred_dat$pred_count <- par_mat[dat$stream,'mu'] * exp(par_mat[dat$stream,'alpha'] * dat$distance)
+
+  pred_dat$ll <- dpois(pred_dat$count,pmax(0,pred_dat$pred_count), log = T)
+
+  #   pred_plot <- pred_dat %>%
+  #   ggplot(aes(distance,count)) +
+  #     geom_point() +
+  #     geom_line(aes(distance, pred_count)) +
+  #     facet_wrap(~stream)
+
+  # Mu prior ----
+
+  mu_expected <- pars[,'mu.f'] + ( (pars[,'mu.l'] - pars[,'mu.f']) * (par_mat[,'stream']-1) ) / (20-1)
+
+  mu_ll <- dnorm(par_mat[,'mu'], mean = mu_expected, sd = pars[,'sigma.mu'], log = T)
+
+  # Alpha prior ----
+
+  alpha_ll <- dnorm(par_mat[,'alpha'],mean = pars[,'alpha.bar'], sd = pars[,'sigma.alpha'], log = T)
+
+  # Hyperprior likelihoods ----
+
+  mu_f_ll <- dunif(pars[,'mu.f'],0,1e10, log = T)
+
+  mu_l_ll <- dunif(pars[,'mu.l'],0,1e10, log = T)
+
+  sigma_mu_ll <- dnorm(pars[,'sigma.mu'], sigma.hyperprior.mu,sigma.hyperprior.sd, log = T)
+
+  sigma_alpha_ll <- dnorm(pars[,'sigma.alpha'], alpha.hyperprior.mu,alpha.hyperprior.sd, log = T)
+
+  ll <- sum(pred_dat$ll) + sum(mu_ll) + sum(alpha_ll) +
+    sigma_mu_ll + sigma_alpha_ll   + mu_f_ll + mu_l_ll
+
+  return(list(loglike = ll,data_loglike = sum(pred_dat$ll), deviance = -2*sum(pred_dat$ll), pred_dat = pred_dat))
+
+}
diff --git a/Fish 558 Homeworks/HWK4/Homework 4 Functions and Data/mozzy_mcmc.R b/Fish 558 Homeworks/HWK4/Homework 4 Functions and Data/mozzy_mcmc.R
@@ -0,0 +1,135 @@
+#' Run MCMC on mosquitoes
+#'
+#' \code{Mozzy_MCMC} runs an MCMC on
+#' mosquito count data
+#' @param par_init initival vector of parameter guesses
+#' @param varnames names of the parameters
+#' @param dat are the raw stream-count-distance data
+#' @param covar is the variance-covariance matrix of pars
+#' @param n_sim number of simulations including burn in
+#' @param n_burn number of simulations to burn in
+#' @param n_thin the number of iterations to thin by
+#' @param vcov_augment factor to multiply vcov by in the burn in period
+#' @param prog_bar T or F to output progress bar
+#' @param seed if you want to set a specific seed for this run
+#' @param targ_accept_rate the target acceptance rate to tune the model
+#' @param jumpyness parameter to modify acceptance ratio
+
+mozzy_mcmc<-function(par_init,parnames,dat,vcov,n_sim=1000,n_burn=0,n_thin=1,
+                     vcov_augment = 1,prog_bar = T,jumpyness = 1, seed = NA,targ_accept_rate = 0.2)
+{
+
+  if (is.na(seed) == F)
+  {
+    set.seed(seed)
+  }
+
+  # Set things up ----
+
+  par_curr <- par_init
+
+  if (class(vcov) == 'data.frame'){ vcov <- as.matrix(vcov) }
+
+  stream_pars = which(is.na(suppressWarnings(as.numeric(gsub('.*\\.', '',parnames)))) == F)
+
+  par_mat <- matrix(NA,nrow = 20,ncol = 3)
+
+  par_mat[,1] <- 1:20
+
+  par_mat[,2:3] <- (par_init[,stream_pars])
+
+  colnames(par_mat) <- c('stream','mu','alpha')
+
+  likelihoods <- mozzy_like(pars = par_curr, par_mat = par_mat, dat = dat,parnames = parnames, stream_pars = stream_pars)
+
+  deviance_curr <- likelihoods$deviance
+
+  ll_curr <- likelihoods$loglike
+
+  outvars <- c(tolower(colnames(par_curr)),'ll','deviance')
+
+  par_out <- matrix(NA, nrow = (n_sim - n_burn), ncol = length(outvars))
+
+  i_count <- 0
+
+  if(prog_bar == T)
+  {
+    pb <- txtProgressBar(min = 1, max = n_sim, style = 3)
+  }
+
+  #   a <- proc.time()
+
+
+  sel_count <- 0
+
+  orig_vcov <- vcov
+
+  vcov_scalar <- vcov_augment
+
+  track_rate <- rep(NA,n_sim)
+
+  ### run MCMC----
+
+  for (i in 1 : n_sim)
+  {
+
+    vcov <- orig_vcov * vcov_scalar
+
+    par_next <- rmvnorm(n = 1, mean= par_curr, sigma =  vcov)
+
+    colnames(par_next) <- parnames
+
+    par_mat[,2:3] <- (par_next[,stream_pars])
+
+    likelihoods <- mozzy_like(pars = par_next, par_mat = par_mat,dat = dat,parnames = parnames,stream_pars = stream_pars)
+
+    deviance_next <- likelihoods$deviance
+
+    ll_next <- likelihoods$loglike
+
+    if(prog_bar == T)
+    {
+
+      setTxtProgressBar(pb, i)
+
+    }
+
+    rand_accept <- log(runif(1,0,jumpyness))
+
+    if (ll_next > ll_curr + rand_accept)
+    {
+      par_curr <- par_next
+      ll_curr <- ll_next
+      deviance_curr <- deviance_next
+      sel_count <- sel_count + 1
+    }
+
+    track_rate[i] <- sel_count/i
+    if (i > n_burn) # store results
+    {
+      i_count <- i_count + 1
+      par_out[i_count,] <- c(as.matrix(par_curr),ll_curr,deviance_curr)
+    } #else{
+    if (i %% round(.2*n_burn,0) == 0 & i <= n_burn) { #Tune variance covariance
+
+      acceptance_rate <- sel_count/i
+
+      perc_off <- (acceptance_rate/targ_accept_rate) * 1.1
+
+      vcov_scalar <- vcov_augment * perc_off
+
+    }
+
+  }
+
+  # proc.time() - a
+  par_out <- as.data.frame(par_out)
+
+  colnames(par_out) <- outvars
+
+  accepted_runs <- data.frame(num_selected = sel_count, num_simulated = n_sim, perc_selected = sel_count/n_sim,
+                              perc_rejected = 1 - sel_count/n_sim, num_kept = n_sim - n_burn)
+
+  return(list(posteriors = par_out[1:i_count,], initial_par = par_init,
+              accepted_runs = accepted_runs))
+}
diff --git a/Fish 558 Homeworks/HWK4/Homework 4 Functions and Data/multi_mozzy_mcmc.R b/Fish 558 Homeworks/HWK4/Homework 4 Functions and Data/multi_mozzy_mcmc.R
@@ -0,0 +1,52 @@
+#' Run multi-chain MCMC on mosquitoes
+#'
+#' \code{multi_mozzy_mcmc} runs potentially multiple chains
+#'  MCMC in parallel
+#' @param par_guess parameter guesses
+#' @param num_starts the number of different starting guesses you wan
+#' @param jitfactor the multiplier for the jitter function
+#' @param numcores number of cores to use in parallel
+#' @param varnames names of the parameters
+#' @param dat are the raw stream-count-distance data
+#' @param covar is the variance-covariance matrix of pars
+#' @param n_sim number of simulations including burn in
+#' @param n_burn number of simulations to burn in
+#' @param n_thin the number of iterations to thin by
+#' @param vcov_augment factor to multiply vcov by in the burn in period
+#' @param prog_bar T or F to output progress bar
+
+multi_mozzy_mcmc<-function(par_guess,num_starts = 1, jitfactor = 20,numcores = 1,parnames,dat,vcov,n_sim=1000,n_burn=0,n_thin=1,
+                           vcov_augment = 1,prog_bar = F,jumpyness = 1)
+{
+
+  registerDoMC(cores=numcores)
+
+    jitter_start <- function(i,pars, jitfactor) #Function to create list of jittered starting guesses
+  {
+
+    newpars <- pars
+
+    newpars[1,] <- jitter(as.numeric(newpars), factor = jitfactor)
+
+    return(newpars)
+  }
+
+  if (num_starts > 1) {
+starting_pars <- lapply(1:num_starts, jitter_start, pars = par_guess,jitfactor = jitfactor) }else{ #Apply jitters, unless you only want to run this once
+  starting_pars <- lapply(1, jitter_start, pars = par_guess, jitfactor = 0)
+}
+
+    #run jags returns list, take a look and see what it actually looks like
+
+  multi_chain <- foreach(i = 1:length(starting_pars)) %dopar% #Run parallel chains
+     mozzy_mcmc(par_init = starting_pars[[i]], parnames = colnames(par_guess), dat = dat,
+                               vcov = vcov, prog_bar = prog_bar, n_sim = n_sim, n_burn = n_burn, n_thin = n_thin,
+                               vcov_augment = vcov_augment,seed = sample(1:1000,1, replace = T))
+
+  # Return posteriors
+  posteriors <- lapply(seq(along = multi_chain), function(i)    data.frame(run = i, multi_chain[[i]]$posteriors))
+
+
+return(list(posteriors = posteriors,multi_chain = multi_chain))
+
+}
diff --git a/Fish 558 Homeworks/HWK4/Homework 4 Functions and Data/run_mozzy_mcmc.R b/Fish 558 Homeworks/HWK4/Homework 4 Functions and Data/run_mozzy_mcmc.R
@@ -0,0 +1,94 @@
+# Wrapper script to run in Terminal
+
+#Initials ----
+
+sims <- 1.5e6
+
+run <- '3.31'
+
+run_folder <- paste('Results/',run,'/',sep = '')
+
+if (dir.exists(run_folder) == F){dir.create(run_folder, recursive = T)}
+
+set.seed(54321)
+
+library(knitr)
+knitr::opts_chunk$set(fig.path='Figs/', echo=FALSE, warning=FALSE, message=FALSE)
+library(gridExtra)
+library(ggplot2)
+library(plyr)
+library(dplyr)
+library(tidyr)
+library(broom)
+library(coda)
+library(ggmcmc)
+# library(LaplacesDemon)
+library(foreach)
+library(scales)
+library(stargazer)
+library(mvtnorm)
+library(doMC)
+library(proftools)
+source('mozzy_likelihood.R')
+source('mozzy_mcmc.R')
+source('multi_mozzy_mcmc.R')
+source('thin_mcmc.R')
+
+# Data ----
+
+dat <- read.csv(file = 'hwk4_data.csv',stringsAsFactors = F) %>%
+  dplyr::rename(stream = Steam) %>%
+  gather('distance','count', grep('X',colnames(.))) %>%
+  mutate(distance = as.numeric(gsub('X','',distance)))
+
+colnames(dat) <- tolower(colnames(dat))
+
+par_guess <- (read.csv(file = 'hwk4_pars.csv', stringsAsFactors = F))
+
+colnames(par_guess) <- tolower(colnames(par_guess))
+
+
+vcov <- read.csv(file = 'hwk4_vcov.csv', stringsAsFactors = F)
+
+rownames(vcov) <- colnames(vcov)
+
+
+dat_plot <- dat %>%
+  ggplot(aes(distance, count, fill = stream, size = count)) +
+  geom_point(shape = 21) +
+  facet_wrap(~stream)
+
+vcov.plot <- vcov %>%
+  mutate(var1 = colnames(.)) %>%
+  gather('var2','covar',which(grepl('var1',colnames(.)) == F)) %>%
+  ggplot(aes(var1,var2,fill = covar)) +
+  geom_tile()
+
+# MCMC ----
+
+
+#
+# optim(.3,tune_scalar,n_sim = n_sim,par_init = as.matrix(par_guess),parnames = colnames(par_guess),dat = dat,vcov = vcov, n_burn=round(.5*n_sim,0)
+#       ,n_thin=1,prog_bar = F,jumpyness = 1, seed = NA,targ_accept_rate = 0.3, lower = .05, upper = 1.5)
+#
+# tune_scalar(.3,n_sim = 1000,par_init = as.matrix(par_guess),parnames = colnames(par_guess),dat = dat,vcov = vcov, n_burn=round(.5*n_sim,0)
+# ,n_thin=1,prog_bar = F,jumpyness = 1, seed = NA,targ_accept_rate = 0.3)
+
+# a <- proc.time()
+
+mcmc_results <- mozzy_mcmc(par_init = as.matrix(par_guess),
+                           parnames = colnames(par_guess), dat = dat,
+                           vcov = vcov, prog_bar = T, n_sim = sims, n_burn = round(.6*sims,0),
+                           n_thin = 1,vcov_augment = (2.4/sqrt(45))^2,jumpyness = 1,targ_accept_rate = .2)
+
+# proc.time() - a
+
+
+mcmc_posteriors <- thin_mcmc(chains = mcmc_results$posteriors, thin_every = sims/1000)
+
+ggmcmc(ggs(mcmc(mcmc_posteriors)), file=paste(run_folder,"mozzy_mcmc_diagnostics.pdf", sep = ''))
+
+
+save(file = paste(run_folder,'mozzy_mcmc.Rdata', sep = ''), mcmc_results)
+
+
diff --git a/Fish 558 Homeworks/HWK4/Homework 4 Functions and Data/run_multi_mcmc.R b/Fish 558 Homeworks/HWK4/Homework 4 Functions and Data/run_multi_mcmc.R
@@ -0,0 +1,72 @@
+# Wrapper script to run in terminal ----
+rm(list = ls())
+
+run <- '3.3'
+
+sims <- 1.5e6
+
+run_folder <- paste('Results/',run,'/',sep = '')
+
+if (dir.exists(run_folder) == F){dir.create(run_folder, recursive = T)}
+
+
+# set.seed(54321)
+library(knitr)
+# knitr::opts_chunk$set(fig.path='Figs/', echo=FALSE, warning=FALSE, message=FALSE)
+library(gridExtra)
+library(ggplot2)
+library(plyr)
+library(dplyr)
+library(tidyr)
+library(broom)
+library(coda)
+library(ggmcmc)
+# library(LaplacesDemon)
+library(foreach)
+library(scales)
+library(stargazer)
+library(mvtnorm)
+library(doMC)
+library(proftools)
+source('mozzy_likelihood.R')
+source('mozzy_mcmc.R')
+source('multi_mozzy_mcmc.R')
+source('thin_mcmc.R')
+
+dat <- read.csv(file = 'hwk4_data.csv',stringsAsFactors = F) %>%
+  dplyr::rename(stream = Steam) %>%
+  gather('distance','count', grep('X',colnames(.))) %>%
+  mutate(distance = as.numeric(gsub('X','',distance)))
+
+colnames(dat) <- tolower(colnames(dat))
+
+par_guess <- (read.csv(file = 'hwk4_pars.csv', stringsAsFactors = F))
+
+colnames(par_guess) <- tolower(colnames(par_guess))
+
+
+vcov <- read.csv(file = 'hwk4_vcov.csv', stringsAsFactors = F)
+
+rownames(vcov) <- colnames(vcov)
+
+
+dat_plot <- dat %>%
+  ggplot(aes(distance, count, fill = stream, size = count)) +
+  geom_point(shape = 21) +
+  facet_wrap(~stream)
+
+vcov.plot <- vcov %>%
+  mutate(var1 = colnames(.)) %>%
+  gather('var2','covar',which(grepl('var1',colnames(.)) == F)) %>%
+  ggplot(aes(var1,var2,fill = covar)) +
+  geom_tile()
+
+multichain_mcmc_results <- multi_mozzy_mcmc(par_guess = as.matrix(par_guess),jitfactor = 5,
+                                            num_starts = 3,numcores = 3, parnames = colnames(par_guess),
+                                            dat = dat, vcov = vcov, prog_bar = F, n_sim = sims,
+                                            n_burn = round(.6*sims,0), n_thin = 1,
+                                            vcov_augment = (2.4/sqrt(45))^2,jumpyness = .001)
+
+print('finished parallel chains')
+
+save(file = paste(run_folder,'multi_chain_mcmc.Rdata',sep = ''), multichain_mcmc_results)
diff --git a/Fish 558 Homeworks/HWK4/Homework 4 Functions and Data/thin_mcmc.R b/Fish 558 Homeworks/HWK4/Homework 4 Functions and Data/thin_mcmc.R
@@ -0,0 +1,7 @@
+thin_mcmc <- function(chains, thin_every =  1)
+{
+
+  thinned_chains <- chains[(1:dim(chains)[1]) %% thin_every == 0,] #only select chains divisible by the thinning interval
+
+  return(thinned_chains)
+}
diff --git a/Fish 558 Homeworks/HWK4/mosquito_likelihood.R b/Fish 558 Homeworks/HWK4/mosquito_likelihood.R
diff --git a/Fish 558 Homeworks/HWK4/mozzy_likelihood.R b/Fish 558 Homeworks/HWK4/mozzy_likelihood.R
@@ -0,0 +1,54 @@
+#' A function to return the likelihood for assignmen
+#'
+#' \code{mosquito.like} returns the negative log likelood of the data
+#' @param pars is a list of parameters for the model
+#' @param dat the stream-distance-count data
+#' @param parnames vector of parameter names
+#' @param stream_pars index of parameters that are stream specific
+#' @param dat data
+#' @param sigma.hyperprior.mu
+
+mozzy_like <- function(pars ,parnames, par_mat,stream_pars, dat , sigma.hyperprior.mu = 1,sigma.hyperprior.sd = 0.1,
+                       alpha.hyperprior.mu = .01, alpha.hyperprior.sd = 0.1)
+{
+
+  # Data likelihood ----
+
+  pred_dat <- dat
+
+  pred_dat$pred_count <- par_mat[dat$stream,'mu'] * exp(par_mat[dat$stream,'alpha'] * dat$distance)
+
+  pred_dat$ll <- dpois(pred_dat$count,pmax(0,pred_dat$pred_count), log = T)
+
+  #   pred_plot <- pred_dat %>%
+  #   ggplot(aes(distance,count)) +
+  #     geom_point() +
+  #     geom_line(aes(distance, pred_count)) +
+  #     facet_wrap(~stream)
+
+  # Mu prior ----
+
+  mu_expected <- pars[,'mu.f'] + ( (pars[,'mu.l'] - pars[,'mu.f']) * (par_mat[,'stream']-1) ) / (20-1)
+
+  mu_ll <- dnorm(par_mat[,'mu'], mean = mu_expected, sd = pars[,'sigma.mu'], log = T)
+
+  # Alpha prior ----
+
+  alpha_ll <- dnorm(par_mat[,'alpha'],mean = pars[,'alpha.bar'], sd = pars[,'sigma.alpha'], log = T)
+
+  # Hyperprior likelihoods ----
+
+  mu_f_ll <- dunif(pars[,'mu.f'],0,1e10, log = T)
+
+  mu_l_ll <- dunif(pars[,'mu.l'],0,1e10, log = T)
+
+  sigma_mu_ll <- dnorm(pars[,'sigma.mu'], sigma.hyperprior.mu,sigma.hyperprior.sd, log = T)
+
+  sigma_alpha_ll <- dnorm(pars[,'sigma.alpha'], alpha.hyperprior.mu,alpha.hyperprior.sd, log = T)
+
+  ll <- sum(pred_dat$ll) + sum(mu_ll) + sum(alpha_ll) +
+    sigma_mu_ll + sigma_alpha_ll   + mu_f_ll + mu_l_ll
+
+  return(list(loglike = ll,data_loglike = sum(pred_dat$ll), deviance = -2*sum(pred_dat$ll), pred_dat = pred_dat))
+
+}
diff --git a/Fish 558 Homeworks/HWK4/mozzy_mcmc.R b/Fish 558 Homeworks/HWK4/mozzy_mcmc.R
@@ -11,16 +11,21 @@
 #' @param n_thin the number of iterations to thin by
 #' @param vcov_augment factor to multiply vcov by in the burn in period
 #' @param prog_bar T or F to output progress bar
+#' @param seed if you want to set a specific seed for this run
+#' @param targ_accept_rate the target acceptance rate to tune the model
+#' @param jumpyness parameter to modify acceptance ratio
 
 mozzy_mcmc<-function(par_init,parnames,dat,vcov,n_sim=1000,n_burn=0,n_thin=1,
-                     vcov_augment = 1,prog_bar = T,jumpyness = 1, seed = NA)
+                     vcov_augment = 1,prog_bar = T,jumpyness = 1, seed = NA,targ_accept_rate = 0.2)
 {
 
   if (is.na(seed) == F)
   {
     set.seed(seed)
   }
 
+  # Set things up ----
+
   par_curr <- par_init
 
   if (class(vcov) == 'data.frame'){ vcov <- as.matrix(vcov) }
@@ -35,9 +40,13 @@ mozzy_mcmc<-function(par_init,parnames,dat,vcov,n_sim=1000,n_burn=0,n_thin=1,
 
   colnames(par_mat) <- c('stream','mu','alpha')
 
-  ll_curr <- mozzy_like(pars = par_curr, par_mat = par_mat, dat = dat,parnames = parnames, stream_pars = stream_pars)$loglike
+  likelihoods <- mozzy_like(pars = par_curr, par_mat = par_mat, dat = dat,parnames = parnames, stream_pars = stream_pars)
+
+  deviance_curr <- likelihoods$deviance
+
+  ll_curr <- likelihoods$loglike
 
-  outvars <- c(tolower(colnames(par_curr)),'ll')
+  outvars <- c(tolower(colnames(par_curr)),'ll','deviance')
 
   par_out <- matrix(NA, nrow = (n_sim - n_burn), ncol = length(outvars))
 
@@ -50,52 +59,69 @@ mozzy_mcmc<-function(par_init,parnames,dat,vcov,n_sim=1000,n_burn=0,n_thin=1,
 
   #   a <- proc.time()
 
-  orig_vcov <- vcov
 
   sel_count <- 0
 
+  orig_vcov <- vcov
+
+  vcov_scalar <- vcov_augment
+
+  track_rate <- rep(NA,n_sim)
+
+  ### run MCMC----
+
   for (i in 1 : n_sim)
   {
 
-    if (i < n_burn) {
-      vcov <- orig_vcov * vcov_augment
-    }else{
-      vcov <- orig_vcov
-    }
+    vcov <- orig_vcov * vcov_scalar
 
     par_next <- rmvnorm(n = 1, mean= par_curr, sigma =  vcov)
 
     colnames(par_next) <- parnames
 
-#     if(par_next[,'mu.f'] < 0 ){browser()}
-
     par_mat[,2:3] <- (par_next[,stream_pars])
 
-    ll_next <- mozzy_like(pars = par_next, par_mat = par_mat,dat = dat,parnames = parnames,stream_pars = stream_pars)$loglike
+    likelihoods <- mozzy_like(pars = par_next, par_mat = par_mat,dat = dat,parnames = parnames,stream_pars = stream_pars)
 
-        if(prog_bar == T)
-        {
+    deviance_next <- likelihoods$deviance
 
-    setTxtProgressBar(pb, i)
+    ll_next <- likelihoods$loglike
 
-        }
+    if(prog_bar == T)
+    {
+
+      setTxtProgressBar(pb, i)
+
+    }
 
     rand_accept <- log(runif(1,0,jumpyness))
 
     if (ll_next > ll_curr + rand_accept)
     {
       par_curr <- par_next
       ll_curr <- ll_next
+      deviance_curr <- deviance_next
       sel_count <- sel_count + 1
     }
 
-    #     if (i > n_burn & i %% n_thin == 0) #%% means divisible by
-    if (i > n_burn) #%% means divisible by
+    track_rate[i] <- sel_count/i
+    if (i > n_burn) # store results
     {
       i_count <- i_count + 1
-      par_out[i_count,] <- c(as.matrix(par_curr),ll_curr)
+      par_out[i_count,] <- c(as.matrix(par_curr),ll_curr,deviance_curr)
+    } #else{
+    if (i %% round(.2*n_burn,0) == 0 & i <= n_burn) { #Tune variance covariance
+
+      acceptance_rate <- sel_count/i
+
+      perc_off <- (acceptance_rate/targ_accept_rate) * 1.1
+
+      vcov_scalar <- vcov_augment * perc_off
+
     }
+
   }
+
   # proc.time() - a
   par_out <- as.data.frame(par_out)
 
@@ -105,5 +131,5 @@ mozzy_mcmc<-function(par_init,parnames,dat,vcov,n_sim=1000,n_burn=0,n_thin=1,
                               perc_rejected = 1 - sel_count/n_sim, num_kept = n_sim - n_burn)
 
   return(list(posteriors = par_out[1:i_count,], initial_par = par_init,
-              accepted_runs = accepted_runs ))
+              accepted_runs = accepted_runs))
 }
diff --git a/Fish 558 Homeworks/HWK4/mozzy_mcmc.Rdata b/Fish 558 Homeworks/HWK4/mozzy_mcmc.Rdata
diff --git a/Fish 558 Homeworks/HWK4/multi_mozzy_mcmc.R b/Fish 558 Homeworks/HWK4/multi_mozzy_mcmc.R
@@ -2,8 +2,7 @@
 #'
 #' \code{multi_mozzy_mcmc} runs potentially multiple chains
 #'  MCMC in parallel
-#' mosquito count data
-#' @param par_guess a dataframe of starting parameter vectors
+#' @param par_guess parameter guesses
 #' @param num_starts the number of different starting guesses you wan
 #' @param jitfactor the multiplier for the jitter function
 #' @param numcores number of cores to use in parallel
@@ -16,13 +15,13 @@
 #' @param vcov_augment factor to multiply vcov by in the burn in period
 #' @param prog_bar T or F to output progress bar
 
-multi_mozzy_mcmc<-function(par_guess,num_starts = 1, jitfactor = 1,numcores = 1,parnames,dat,vcov,n_sim=1000,n_burn=0,n_thin=1,
+multi_mozzy_mcmc<-function(par_guess,num_starts = 1, jitfactor = 20,numcores = 1,parnames,dat,vcov,n_sim=1000,n_burn=0,n_thin=1,
                            vcov_augment = 1,prog_bar = F,jumpyness = 1)
 {
 
   registerDoMC(cores=numcores)
 
-    jitter_start <- function(i,pars, jitfactor)
+    jitter_start <- function(i,pars, jitfactor) #Function to create list of jittered starting guesses
   {
 
     newpars <- pars
@@ -33,19 +32,20 @@ multi_mozzy_mcmc<-function(par_guess,num_starts = 1, jitfactor = 1,numcores = 1,
   }
 
   if (num_starts > 1) {
-starting_pars <- lapply(1:num_starts, jitter_start, pars = par_guess,jitfactor = jitfactor) }else{
+starting_pars <- lapply(1:num_starts, jitter_start, pars = par_guess,jitfactor = jitfactor) }else{ #Apply jitters, unless you only want to run this once
   starting_pars <- lapply(1, jitter_start, pars = par_guess, jitfactor = 0)
 }
 
+    #run jags returns list, take a look and see what it actually looks like
 
-  multi_chain <- foreach(i = 1:length(starting_pars)) %dopar%
+  multi_chain <- foreach(i = 1:length(starting_pars)) %dopar% #Run parallel chains
      mozzy_mcmc(par_init = starting_pars[[i]], parnames = colnames(par_guess), dat = dat,
                                vcov = vcov, prog_bar = prog_bar, n_sim = n_sim, n_burn = n_burn, n_thin = n_thin,
-                               vcov_augment = vcov_augment,seed = runif(1,1,1000))
+                               vcov_augment = vcov_augment,seed = sample(1:1000,1, replace = T))
 
+  # Return posteriors
+  posteriors <- lapply(seq(along = multi_chain), function(i)    data.frame(run = i, multi_chain[[i]]$posteriors))
 
-  posteriors <- lapply(seq(along = multi_chain), function(i)    data.frame(run = i, multi_chain[[i]]$posteriors)) %>%
-    ldply()
 
 return(list(posteriors = posteriors,multi_chain = multi_chain))
 

diff --git a/Fish 558 Homeworks/HWK4/posterior_predictive_mozzies.R b/Fish 558 Homeworks/HWK4/posterior_predictive_mozzies.R
@@ -0,0 +1,12 @@
+posterior_predictive_mozzies <- function()
+{
+
+  #1. take the posteriors from the mcmc
+  #2. generate Mul, MuF, alpha, etc. from their appropriate distributions, with para
+  #meters drawn from posterior
+  #3. Generate u_i from the normal distribution with the above parameters
+  #4. predict each count along the river
+
+
+
+}
diff --git a/Fish 558 Homeworks/HWK4/run_mozzy_mcmc.R b/Fish 558 Homeworks/HWK4/run_mozzy_mcmc.R
@@ -0,0 +1,94 @@
+# Wrapper script to run in Terminal
+
+#Initials ----
+
+sims <- 1.5e6
+
+run <- '3.31'
+
+run_folder <- paste('Results/',run,'/',sep = '')
+
+if (dir.exists(run_folder) == F){dir.create(run_folder, recursive = T)}
+
+set.seed(54321)
+
+library(knitr)
+knitr::opts_chunk$set(fig.path='Figs/', echo=FALSE, warning=FALSE, message=FALSE)
+library(gridExtra)
+library(ggplot2)
+library(plyr)
+library(dplyr)
+library(tidyr)
+library(broom)
+library(coda)
+library(ggmcmc)
+# library(LaplacesDemon)
+library(foreach)
+library(scales)
+library(stargazer)
+library(mvtnorm)
+library(doMC)
+library(proftools)
+source('mozzy_likelihood.R')
+source('mozzy_mcmc.R')
+source('multi_mozzy_mcmc.R')
+source('thin_mcmc.R')
+
+# Data ----
+
+dat <- read.csv(file = 'hwk4_data.csv',stringsAsFactors = F) %>%
+  dplyr::rename(stream = Steam) %>%
+  gather('distance','count', grep('X',colnames(.))) %>%
+  mutate(distance = as.numeric(gsub('X','',distance)))
+
+colnames(dat) <- tolower(colnames(dat))
+
+par_guess <- (read.csv(file = 'hwk4_pars.csv', stringsAsFactors = F))
+
+colnames(par_guess) <- tolower(colnames(par_guess))
+
+
+vcov <- read.csv(file = 'hwk4_vcov.csv', stringsAsFactors = F)
+
+rownames(vcov) <- colnames(vcov)
+
+
+dat_plot <- dat %>%
+  ggplot(aes(distance, count, fill = stream, size = count)) +
+  geom_point(shape = 21) +
+  facet_wrap(~stream)
+
+vcov.plot <- vcov %>%
+  mutate(var1 = colnames(.)) %>%
+  gather('var2','covar',which(grepl('var1',colnames(.)) == F)) %>%
+  ggplot(aes(var1,var2,fill = covar)) +
+  geom_tile()
+
+# MCMC ----
+
+
+#
+# optim(.3,tune_scalar,n_sim = n_sim,par_init = as.matrix(par_guess),parnames = colnames(par_guess),dat = dat,vcov = vcov, n_burn=round(.5*n_sim,0)
+#       ,n_thin=1,prog_bar = F,jumpyness = 1, seed = NA,targ_accept_rate = 0.3, lower = .05, upper = 1.5)
+#
+# tune_scalar(.3,n_sim = 1000,par_init = as.matrix(par_guess),parnames = colnames(par_guess),dat = dat,vcov = vcov, n_burn=round(.5*n_sim,0)
+# ,n_thin=1,prog_bar = F,jumpyness = 1, seed = NA,targ_accept_rate = 0.3)
+
+# a <- proc.time()
+
+mcmc_results <- mozzy_mcmc(par_init = as.matrix(par_guess),
+                           parnames = colnames(par_guess), dat = dat,
+                           vcov = vcov, prog_bar = T, n_sim = sims, n_burn = round(.6*sims,0),
+                           n_thin = 1,vcov_augment = (2.4/sqrt(45))^2,jumpyness = 1,targ_accept_rate = .2)
+
+# proc.time() - a
+
+
+mcmc_posteriors <- thin_mcmc(chains = mcmc_results$posteriors, thin_every = sims/1000)
+
+ggmcmc(ggs(mcmc(mcmc_posteriors)), file=paste(run_folder,"mozzy_mcmc_diagnostics.pdf", sep = ''))
+
+
+save(file = paste(run_folder,'mozzy_mcmc.Rdata', sep = ''), mcmc_results)
+
+
diff --git a/Fish 558 Homeworks/HWK4/run_multi_mcmc.R b/Fish 558 Homeworks/HWK4/run_multi_mcmc.R
@@ -1,4 +1,15 @@
+# Wrapper script to run in terminal ----
 rm(list = ls())
+
+run <- '3.3'
+
+sims <- 1.5e6
+
+run_folder <- paste('Results/',run,'/',sep = '')
+
+if (dir.exists(run_folder) == F){dir.create(run_folder, recursive = T)}
+
+
 # set.seed(54321)
 library(knitr)
 # knitr::opts_chunk$set(fig.path='Figs/', echo=FALSE, warning=FALSE, message=FALSE)
@@ -17,7 +28,7 @@ library(stargazer)
 library(mvtnorm)
 library(doMC)
 library(proftools)
-source('mosquito_likelihood.R')
+source('mozzy_likelihood.R')
 source('mozzy_mcmc.R')
 source('multi_mozzy_mcmc.R')
 source('thin_mcmc.R')
@@ -50,11 +61,12 @@ vcov.plot <- vcov %>%
   ggplot(aes(var1,var2,fill = covar)) +
   geom_tile()
 
-sims <- 1e6
-
-
-mcmc_results <- multi_mozzy_mcmc(par_guess = as.matrix(par_guess),num_starts = 3,numcores = 3, parnames = colnames(par_guess), dat = dat, vcov = vcov, prog_bar = F, n_sim = sims, n_burn = round(.5*sims,0), n_thin = 1,vcov_augment = .2,jumpyness = .001)
+multichain_mcmc_results <- multi_mozzy_mcmc(par_guess = as.matrix(par_guess),jitfactor = 5,
+                                            num_starts = 3,numcores = 3, parnames = colnames(par_guess),
+                                            dat = dat, vcov = vcov, prog_bar = F, n_sim = sims,
+                                            n_burn = round(.6*sims,0), n_thin = 1,
+                                            vcov_augment = (2.4/sqrt(45))^2,jumpyness = .001)
 
 print('finished parallel chains')
 
-save(file = 'multi_chain_mcmc.Rdata', mcmc_results)
+save(file = paste(run_folder,'multi_chain_mcmc.Rdata',sep = ''), multichain_mcmc_results)
diff --git a/Fish 558 Homeworks/HWK4/tune_scalar.R b/Fish 558 Homeworks/HWK4/tune_scalar.R
@@ -0,0 +1,12 @@
+tune_scalar <- function(scalar,n_sim = 10000,par_init,parnames,dat,vcov,n_burn=0,n_thin=1,
+                     vcov_augment = 1,prog_bar = T,jumpyness = 1, seed = NA,targ_accept_rate = 0.3)
+{
+
+  mcmc_results <- mozzy_mcmc(par_init = par_init,
+                             parnames = parnames, dat = dat,
+                             vcov = vcov, prog_bar = prog_bar, n_sim = n_sim, n_burn = round(.5*n_sim,0),
+                             n_thin = 1,vcov_augment = scalar,jumpyness = .001,targ_accept_rate = targ_accept_rate)
+
+  return((mcmc_results$accepted_runs$perc_selected - targ_accept_rate)^2)
+
+}
diff --git a/Fish 558 Homeworks/HWK5/HWK5.Rmd b/Fish 558 Homeworks/HWK5/HWK5.Rmd
diff --git a/Fish 558 Homeworks/HWK5/HWK5_Shiny.Rmd b/Fish 558 Homeworks/HWK5/HWK5_Shiny.Rmd
@@ -0,0 +1,44 @@
+---
+title: "FISH-558-Shiny-Homework-5"
+author: "Dan Ovando"
+date: "November 13, 2015"
+output: html_document
+runtime: shiny
+---
+
+
+The goal of this assignment is to provide and MSE/decision analysis around the management of a mythical endangered species. There is uncertainty around parameters and a variety of potential management actions that need considering. 
+
+```{r load parameters}
+
+# n_at_age <- matrix(NA,nrow = sim_years, ncol = ages)
+
+```
+
+## Task A
+
+Find the relationship between the mean and CB of a beta distribtuion and the shape1 and shape 2 parameters for rbeta. The goal here is to take the mean and CV parameters for the beta distribution provided by Andre and convert them into R usable form. 
+
+The beta distribution is a continuous probability distribution defined on the interval [0,1], that can take a wide array of shapes 
+
+```{r, echo=FALSE}
+inputPanel(
+  selectInput("n_breaks", label = "Number of bins:",
+              choices = c(10, 20, 35, 50), selected = 20),
+  
+  sliderInput("bw_adjust", label = "Bandwidth adjustment:",
+              min = 0.2, max = 2, value = 1, step = 0.2)
+)
+
+renderPlot({
+  hist(faithful$eruptions, probability = TRUE, breaks = as.numeric(input$n_breaks),
+       xlab = "Duration (minutes)", main = "Geyser eruption duration")
+  
+  dens <- density(faithful$eruptions, adjust = input$bw_adjust)
+  lines(dens, col = "blue")
+})
+```
+
+
+
+
diff --git a/Fish 558 Homeworks/HWK5/HWK5_cache/html/__packages b/Fish 558 Homeworks/HWK5/HWK5_cache/html/__packages
@@ -0,0 +1,17 @@
+base
+shiny
+knitr
+gridExtra
+ggplot2
+plyr
+dplyr
+tidyr
+broom
+foreach
+scales
+stargazer
+mvtnorm
+proftools
+iterators
+doMC
+DT
diff --git a/Fish 558 Homeworks/HWK5/moments_to_shape.R b/Fish 558 Homeworks/HWK5/moments_to_shape.R
@@ -0,0 +1,21 @@
+#' Convert moment to shape parameters
+#'
+#' \code{moments_to_shape} takes the mean and cv
+#' of a beta distribution and converts it to the shape
+#' parameters
+#' @param mu the mean of the beta distribution
+#' @param cv the cv of the beta distribution
+
+moments_to_shape <- function(mu,cv)
+{
+
+  var <- (mu*cv)^2 #calculate variance
+
+  a <- ( (1-mu)/var - 1/mu )*mu^2
+
+  b <- a*(1/mu - 1)
+
+  if (mu >=1){warning('mu is greater than or equal to 1, reduce it!')}
+
+  return(list(shape1 = a, shape2 = b))
+}
diff --git a/Fish 558 Homeworks/HWK5/sim_problembeast.R b/Fish 558 Homeworks/HWK5/sim_problembeast.R
@@ -0,0 +1,143 @@
+#' runs problembeast model
+#'
+#'\code{sim_problembeast} simulates populations of
+#'the endangered problembeast forward under
+#'different states of nature and hunting
+#'@param pb a list of problembeast population
+#'parameters
+#'@param sim_years number of years to run
+#'@param hunted_patch the patch in which hunting
+#'occurs
+#'@param s_j juvenile survival
+#'@param _value the price per hunted animal
+#'at various ages
+#'@param results blank matrix to store results
+#'@param recruits index of recruits
+#'@param juveniles index of juveniles
+#'@param adults index of all adults
+#'@param mature index of mature adults
+#'@param plus index of the plus group
+#'
+sim_problembeast <- function(pb,sim_years = 101,hunted_patch = 0,s_j = .75,
+                             recruit_value = 1, juv_value = 5, adult_value = 30,
+                             results, recruits,juveniles,adults,mature,plus)
+{
+
+pb$s_j <- s_j
+
+# resultnames <- c('year','adults','hunted','hunt_value','hunted_patch','s_j')
+#
+# results <- as.data.frame(matrix(0,nrow = sim_years+1,ncol = length(resultnames)))
+#
+# colnames(results) <- resultnames
+
+results$year <- 1:(sim_years+1)
+
+results$hunted_patch <- hunted_patch
+
+results$s_j <- s_j
+
+# recruits <- which(colnames(pb$n_at_age) %in% 'age.0')
+#
+# juveniles <- which(colnames(pb$n_at_age) %in% paste('age',1:8,sep = '.'))
+#
+# adults <- which(colnames(pb$n_at_age) %in% paste('age',9:14,sep = '.'))
+#
+# mature <- which(colnames(pb$n_at_age) %in% paste('age',9:15,sep = '.'))
+#
+# plus <- which(colnames(pb$n_at_age) == 'age.15')
+
+phi <- matrix(1, nrow = 3, ncol = sim_years)
+
+if (hunted_patch>0){
+
+  shape1 <- eval(parse(text = paste('pb$phi_',hunted_patch,'_shape',sep = '') ))$shape1
+
+  shape2 <- eval(parse(text = paste('pb$phi_',hunted_patch,'_shape',sep = '') ))$shape2
+
+  phi[hunted_patch,] <-  rbeta(sim_years,shape1 = shape1, shape2 = shape2)
+
+}
+
+n_adults <- sum(pb$n_at_age[1, mature])
+
+results$adults[1] <- n_adults
+
+pb$n_at_age <- as.matrix(pb$n_at_age)
+
+  for (y in 1:(sim_years))
+  {
+
+    hunted <- 0
+
+    juv_survive <- rbinom(length(juveniles),(pb$n_at_age[y,juveniles - 1]),pb$s_j)
+
+    juv_survive_hunting <- rbinom(length(juveniles),(juv_survive),phi[2,y])
+
+    juv_hunted <- juv_survive - juv_survive_hunting
+
+    pb$n_at_age[y + 1,juveniles] <- juv_survive_hunting
+
+    adults_survive <- rbinom(length(adults),(pb$n_at_age[y,adults - 1]),pb$s_a)
+
+    adults_survive_hunting <- rbinom(length(adults), (adults_survive),phi[3,y])
+
+    adults_hunted <- adults_survive - adults_survive_hunting
+
+    pb$n_at_age[y + 1, adults] <- adults_survive_hunting
+
+    enter_plus_survive <- rbinom(length(plus),(pb$n_at_age[y,plus - 1]),pb$s_a)
+
+    enter_plus_survive_hunting <- rbinom(length(plus), (enter_plus_survive),phi[3,y])
+
+    enter_plus_hunted <- enter_plus_survive - enter_plus_survive_hunting
+
+    enter_plus <- enter_plus_survive_hunting
+
+    stay_plus_survive <- rbinom(length(plus),(pb$n_at_age[y,plus]),pb$s_a)
+
+    stay_plus_survive_hunted <- rbinom(length(plus), (stay_plus_survive),phi[3,y])
+
+    stay_plus_hunted <- stay_plus_survive - stay_plus_survive_hunted
+
+    stay_plus <- stay_plus_survive_hunted
+
+    hunted_in_two <- sum(juv_hunted)
+
+    hunted_in_three <- sum(adults_hunted) + sum(enter_plus_hunted) + sum(stay_plus_hunted)
+
+    hunted <-  hunted_in_two + hunted_in_three
+
+    hunted_value <- hunted_in_two * juv_value + hunted_in_three * adult_value
+
+    pb$n_at_age[y + 1, plus] <- enter_plus + stay_plus
+
+    n_adults <- sum(pb$n_at_age[y + 1, mature])
+
+    recruits_survive <- rbinom(length(recruits), n_adults, pb$b)
+
+    recruits_survive_hunting <- rbinom(length(recruits), (recruits_survive),phi[1,y])
+
+    recruits_hunted <- recruits_survive - recruits_survive_hunting
+
+    if (hunted_patch >1){
+    results$hunted[y] <- hunted
+    results$hunt_value[y] <- hunted_value
+    }
+    if(hunted_patch == 1)
+    {
+      results$hunted[y + 1] <- recruits_hunted
+      results$hunt_value[y+1] <- recruits_hunted * recruit_value
+    }
+
+    pb$n_at_age[y + 1,recruits] <- recruits_survive_hunting
+
+    results$adults[y+1] <- n_adults
+
+  }
+
+results <- results[1:sim_years,]
+
+return(results)
+
+}
diff --git a/Fish 558 Homeworks/Homework 2/HWK2.doc b/Fish 558 Homeworks/Homework 2/HWK2.doc
diff --git a/Fish 558 Workshop/.DS_Store b/Fish 558 Workshop/.DS_Store