mbjoseph

---
title: "Untitled"
author: "Max Joseph"
date: "2/6/2019"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}

mbjoseph

mbjoseph

# Simulation script for factor analysis ala Leung & Drton (2016) ----------
library(rstan)
library(bayesplot)

m <- 5 # dimension of observed data (e.g., # traits)
k <- 2 # number of latent factors
n <- 100 # number of sample units (e.g., # species)

# residual variance matrix (is diagonal)
Omega <- diag(.3 + abs(rnorm(m, sd = .3)))

mbjoseph

# Script to get MCD14ML MODIS data
library(RCurl)
library(dplyr)

# fetch_MCD14ML() downloads all of the zipped data files for the
# MODIS MCD14ML data product
#
# args:
#   dir: the name of a directory to store the zip files (string)
#   overwrite: should existing data be overrwitten? (logical)

mbjoseph

library(lme4)

# adapted from ?predict.merMod
gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 |herd), 
             cbpp, 
             binomial)

# add some new levels to herd, which is a random effect
newdata <- with(cbpp, 
                expand.grid(period=unique(period), 

mbjoseph

library(tidyverse)
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

# Simulate true occupancy states ------------------------------------------
# define a design matrix for site-level occupancy
n_site <- 200
m_psi <- 3    # m_psi is the number of columns in the site level design matrix
X_psi <- matrix(c(rep(1, n_site), rnorm((m_psi - 1) * n_site)), 

mbjoseph

library(tidyverse)
library(reshape2)
library(pbapply)

n_samples <- 100
n_wavelength <- 400

wavelengths <- matrix(nrow = n_samples, ncol = n_wavelength)
reflectances <- matrix(nrow = n_samples, ncol = n_wavelength)
for (i in 1:n_samples) {

mbjoseph

mbjoseph

# Distance sampling model with data augmentation. 
# Based on section 8.3.1 of the Applied Hierarchical Modeling book by Royle and Kery

library(rstan)

B <- 50

# note that I'm dividing by 10 and adding 1e-6 to put this on a manageable
# scale and prevent values == 0
d_obs <- c(71.93, 26.05, 58.47, 92.35, 163.83, 84.52, 163.83, 157.33,

mbjoseph

vioplot2 <- function (x, ..., range = 1.5, h = NULL, ylim = NULL, names = NULL, 
                      horizontal = FALSE, col = "magenta", border = "black", lty = 1, 
                      lwd = 1, rectCol = "black", colMed = "white", pchMed = 19, 
                      at, add = FALSE, wex = 1, drawRect = TRUE, side="both") 
{
  datas <- list(x, ...)
  n <- length(datas)
  if (missing(at)) 
    at <- 1:n
  upper <- vector(mode = "numeric", length = n)