trying methods for looping over replicates efficiently

inst/examples/myers-exploration.Rmd

@@ -0,0 +1,67 @@
+```{r set-options, echo = FALSE, cache = FALSE, external = TRUE, include = FALSE}
+opts_chunk$set(external = TRUE, cache = FALSE, cache.path = "myers-cache/", warning=FALSE)
+read_chunk('gaussian-process-control.R')
+library(knitcitations)
+```
+
+
+```{r libraries, include=FALSE}
+```
+```{r graphing-options, include=FALSE}
+```
+
+
+Fixed priors on hyperparameters, fixed model type.
+
+```{r gp-priors}
+```
+
+```{r sdp-pars-fixed}
+profit = function(x,h) pmin(x, h)
+delta <- 0.01
+OptTime = 20
+reward = 0
+xT <- 0
+z_g = function() rlnorm(1, 0, sigma_g)
+z_m = function() 1+(2*runif(1, 0,  1)-1) * sigma_m
+```
+
+
+```{r Myer-explore}
+f <- Myer_harvest
+pars <- c(1, 2, 4.5) 
+p <- pars # shorthand 
+K <- p[1] * p[3] / 2 + sqrt( (p[1] * p[3]) ^ 2 - 4 * p[3] ) / 2
+allee <- p[1] * p[3] / 2 - sqrt( (p[1] * p[3]) ^ 2 - 4 * p[3] ) / 2 # allee threshold
+e_star <- (p[1] * sqrt(p[3]) - 2) / 2 ## Bifurcation point 
+```
+
+
+```{r sdp-pars-explore}
+sigma_g <- 0.05
+sigma_m <- 0.2
+x_grid <- seq(0, 1.5 * K, length=101)
+h_grid <- x_grid
+```
+
+With parameters `r p`. 
+
+
+```{r }
+yields <- 
+  lapply(1:3, function(j){
+  source("single-gp-fit.R")
+  yield
+})
+```
+
+```{r}
+yields <- melt(yields, id=c("method", "V1", "sd"))
+yields         
+````
+
+
+
+```{r echo=FALSE, results="asis"}
+bibliography("html")
+```

inst/examples/single-gp-fit.R

@@ -0,0 +1,143 @@
+x_0_observed <- allee + x_grid[5]
+xT <- 0
+seed <- round(runif(1) * 1e6)
+seed
+set.seed(seed)
+
+## @knitr sim-obs
+Tobs <- 50
+x <- numeric(Tobs)
+x[1] <- x_0_observed
+for(t in 1:(Tobs-1))
+  x[t+1] = z_g() * f(x[t], h=0, p=p)
+plot(x)
+
+## @knitr lag-data
+obs <- data.frame(x=c(0,x[1:(Tobs-1)]),y=c(0,x[2:Tobs]))
+
+## @knitr par-est
+estf <- function(p){
+  mu <- log(obs$x) + p["r"]*(1-obs$x/p["K"])
+  -sum(dlnorm(obs$y, mu, p["s"]), log=TRUE)
+}
+o <- optim(par = c(r=1,K=mean(x),s=1), estf, method="L", lower=c(1e-3,1e-3,1e-3))
+f_alt <- Ricker
+p_alt <- c(o$par['r'], o$par['K'])
+sigma_g_alt <- o$par['s']
+
+
+## @knitr gp-fit
+gp <- bgp(X=obs$x, XX=x_grid, Z=obs$y, verb=0,
+          meanfn="constant", bprior="b0", BTE=c(2000,16000,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")
+
+
+## @knitr gp-data
+V <- gp$ZZ.ks2
+Ef = gp$ZZ.km
+tgp_dat <- data.frame(x   = gp$XX[[1]], 
+                      y   = gp$ZZ.km, 
+                      ymin = gp$ZZ.km - 1.96 * sqrt(gp$ZZ.ks2), 
+                      ymax = gp$ZZ.km + 1.96 * sqrt(gp$ZZ.ks2))
+
+
+## @knitr gp-plot
+true <- sapply(x_grid, f, 0, p)
+est <- sapply(x_grid, f_alt, 0, p_alt)
+models <- data.frame(x=x_grid, GP=tgp_dat$y, Parametric=est, True=true)
+models <- melt(models, id="x")
+names(models) <- c("x", "method", "value")
+ggplot(tgp_dat)  + geom_ribbon(aes(x,y,ymin=ymin,ymax=ymax), fill="gray80") +
+  geom_line(data=models, aes(x, value, col=method), lwd=2, alpha=0.8) + 
+  geom_point(data=obs, aes(x,y), alpha=0.8) + 
+  xlab(expression(X[t])) + ylab(expression(X[t+1])) +
+  scale_colour_manual(values=cbPalette)
+
+
+
+## @knitr gp-posteriors
+hyperparameters <- c("index", "s2", "tau2", "beta0", "nug", "d", "ldetK")
+posteriors <- melt(gp$trace$XX[[1]][,hyperparameters], id="index")
+priors <- list(s2 = s2_prior, tau2 = tau2_prior, beta0 = dnorm, nug = nug_prior, d = d_prior, ldetK = function(x) 0)
+prior_curves <- ddply(posteriors, "variable", function(dd){
+  grid <- seq(min(dd$value), max(dd$value), length = 100)
+  data.frame(value = grid, density = priors[[dd$variable[1]]](grid))
+})
+ggplot(posteriors) + 
+  geom_histogram(aes(x=value, y=..density..), lwd=2) +
+  geom_line(data=prior_curves, aes(x=value, y=density), col="red", lwd=2) +
+  facet_wrap(~ variable, scale="free")
+ggplot(prior_curves) + 
+  geom_line(aes(x=value, y=density), col="red", lwd=2) +
+  facet_wrap(~ variable, scale="free")
+
+## @knitr gp-opt
+matrices_gp <- gp_transition_matrix(Ef, V, x_grid, h_grid)
+opt_gp <- find_dp_optim(matrices_gp, x_grid, h_grid, OptTime, xT, profit, delta, reward=reward)
+
+## @knitr true-opt
+matrices_true <- f_transition_matrix(f, p, x_grid, h_grid, sigma_g)
+opt_true <- find_dp_optim(matrices_true, x_grid, h_grid, OptTime, xT, profit, delta=delta, reward = reward)
+
+## @knitr est-opt
+matrices_estimated <- f_transition_matrix(f_alt, p_alt, x_grid, h_grid, sigma_g_alt)
+opt_estimated <- find_dp_optim(matrices_estimated, x_grid, h_grid, OptTime, xT, profit, delta=delta, reward = reward)
+
+## @knitr policy_plot
+policies <- melt(data.frame(stock=x_grid, 
+                            GP = x_grid[opt_gp$D[,1]], 
+                            Parametric = x_grid[opt_estimated$D[,1]],
+                            True = x_grid[opt_true$D[,1]]),
+                 id="stock")
+names(policies) <- c("stock", "method", "value")
+policy_plot <- ggplot(policies, aes(stock, stock - value, color=method)) +
+  geom_line(lwd=2, alpha=0.8) + 
+  xlab("stock size") + ylab("escapement")  +
+  scale_colour_manual(values=cbPalette)
+policy_plot
+
+
+
+## @knitr stationary_policy_only
+m <- sapply(1:OptTime, function(i) opt_gp$D[,1])
+opt_gp$D <- m
+mm <- sapply(1:OptTime, function(i) opt_true$D[,1])
+opt_true$D <- mm
+mmm <- sapply(1:OptTime, function(i) opt_estimated$D[,1])
+opt_estimated$D <- mmm
+
+## @knitr simulate
+set.seed(1)
+sim_gp <- lapply(1:100, function(i) ForwardSimulate(f, p, x_grid, h_grid, K, opt_gp$D, z_g, profit=profit))
+set.seed(1)
+sim_true <- lapply(1:100, function(i) ForwardSimulate(f, p, x_grid, h_grid, K, opt_true$D, z_g, profit=profit))
+set.seed(1)
+sim_est <- lapply(1:100, function(i) ForwardSimulate(f, p, x_grid, h_grid, K, opt_estimated$D, z_g, profit=profit))
+
+
+## @knitr tidy
+dat <- list(GP = sim_gp, Parametric = sim_est, True = sim_true)
+dat <- melt(dat, id=names(dat[[1]][[1]]))
+dt <- data.table(dat)
+setnames(dt, c("L1", "L2"), c("method", "reps")) 
+
+## @knitr sim-fish
+ggplot(dt) + 
+  geom_line(aes(time, fishstock, group=interaction(reps,method), color=method), alpha=.1) +
+  scale_colour_manual(values=cbPalette, guide = guide_legend(override.aes = list(alpha = 1)))
+
+## @knitr sim-harvest
+ggplot(dt) +
+  geom_line(aes(time, harvest, group=interaction(reps,method), color=method), alpha=.1) +
+  scale_colour_manual(values=cbPalette, guide = guide_legend(override.aes = list(alpha = 1)))
+
+
+## @knitr costs
+profits <- dt[, sum(profit), by = c("reps", "method")]
+means <- profits[, mean(V1), by = method]
+sds <- profits[, sd(V1), by = method]
+yield <- cbind(means, sd = sds$V1)
+yield
+