opts_chunk$set(external = TRUE, cache = FALSE, cache.path = "allee/", warning=FALSE) library(knitcitations)
require(pdgControl) require(nonparametricbayes) require(reshape2) require(ggplot2) require(data.table) require(tgp) library(kernlab) require(MCMCpack) require(plyr)
opts_knit$set(upload.fun = socialR::flickr.url) #opts_chunk$set(dev = 'Cairo_pdf', dev.args=list("")) opts_chunk$set(dev="png", dev.args=list(bg="transparent")) opts_chunk$set(comment=NA, tidy=FALSE) theme_set(theme_bw(base_size=16)) theme_update(panel.background = element_rect(fill = "transparent",colour = NA), plot.background = element_rect(fill = "transparent",colour = NA)) cbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
Fixed priors on hyperparameters, fixed model type.
#inv gamma has mean b / (a - 1) (assuming a>1) and variance b ^ 2 / ((a - 2) * (a - 1) ^ 2) (assuming a>2) s2.p <- c(5,5) tau2.p <- c(5,1) d.p = c(10, 1/0.1, 10, 1/0.1) nug.p = c(10, 1/0.1, 10, 1/0.1) # gamma mean s2_prior <- function(x) dinvgamma(x, s2.p[1], s2.p[2]) tau2_prior <- function(x) dinvgamma(x, tau2.p[1], tau2.p[2]) d_prior <- function(x) dgamma(x, d.p[1], scale = d.p[2]) + dgamma(x, d.p[3], scale = d.p[4]) nug_prior <- function(x) dgamma(x, nug.p[1], scale = nug.p[2]) + dgamma(x, nug.p[3], scale = nug.p[4]) beta0_prior <- function(x, tau) dnorm(x, 0, tau) beta = c(0) priors <- list(s2 = s2_prior, tau2 = tau2_prior, beta0 = dnorm, nug = nug_prior, d = d_prior, ldetK = function(x) 0)
profit = function(x,h) pmin(x, h) delta <- 0.01 OptTime = 20 # stationarity with unstable models is tricky thing reward = 0 xT <- 0
f <- RickerAllee # c(5, 10, 5) is 2-cycle, c(5.5, 10, 5) is 6 cycle, 5.3 is about 4 p <- c(2, 10, 5) K <- 10 allee <- 5
sigma_g <- 0.05 sigma_m <- 0.0 z_g = function() rlnorm(1, 0, sigma_g) z_m = function() 1+(2*runif(1, 0, 1)-1) * sigma_m x_grid <- seq(0, 1.5 * K, length=101) h_grid <- x_grid
With parameters r p
.
seed_i <- 1 Xo <- K # observations start from x0 <- Xo # simulation under policy starts from
obs <- sim_obs(Xo, z_g, f, p, Tobs=35, nz=15, harvest = sort(rep(seq(0, .5, length=7), 5)), seed = seed_i)
alt <- par_est(obs, init = c(r=p[1], K=mean(obs$x), s=sigma_g)) est <- par_est_allee(obs, f, p, init = c(2, mean(obs$x), 2, s = sigma_g))
Which estimates a Ricker model with $r =$ r alt$p[1]
, $K =$ r alt$p[2]
, and the Allen allee model with $r =$ r est$p[1]
, $K =$ r est$p[2]
and $C =$ r est$p[3]
.
gp <- bgp(X=obs$x, XX=x_grid, Z=obs$y, verb=0, meanfn="constant", bprior="b0", BTE=c(10,1600,2), m0r1=FALSE, corr="exp", trace=TRUE, beta = beta, s2.p = s2.p, d.p = d.p, nug.p = nug.p, tau2.p = tau2.p, s2.lam = "fixed", d.lam = "fixed", nug.lam = "fixed", tau2.lam = "fixed") gp_plot(gp, f, p, est$f, est$p, alt$f, alt$p, x_grid, obs, seed_i)
posteriors_plot(gp, priors) # needs trace=TRUE!
OPT <- optimal_policy(gp, f, est$f, alt$f, p, est$p, alt$p, x_grid, h_grid, sigma_g, sigma_g, sigma_g, # est$sigma_g, alt$sigma_g, but those ests are poor delta, xT, profit, reward, OptTime) plot_policies(x_grid, OPT$gp_D, OPT$est_D, OPT$true_D, OPT$alt_D)
dt <- simulate_opt(OPT, f, p, x_grid, h_grid, x0, z_g, profit, OptTime=OptTime) sim_plots(dt, seed=seed_i) profits_stats(dt)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.