Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

added value of information example

  • Loading branch information...
commit e5d913b9e06822e83702f0f41cd67275d3a1e533 1 parent 80c65d5
@cboettig authored
View
14 R/ForwardSimulate.R
@@ -15,11 +15,13 @@
#' assessment of stock size
#' @param z_i a function which returns a random number from a distribution
#' for the implementation uncertainty in quotas
+#' @param profit function, if known
#' @return a data frame with the time, fishstock, harvested amount,
#' and what the escapement ("unharvested").
#' @export
ForwardSimulate <- function(f, pars, x_grid, h_grid, x0, D, z_g,
- z_m=function(x) 1, z_i = function(x) 1){
+ z_m=function(x) 1, z_i = function(x) 1,
+ profit=NULL){
# initialize variables with initial conditions
OptTime <- dim(D)[2] # Stopping time
x_h <- numeric(OptTime) # population dynamics with harvest
@@ -28,7 +30,8 @@ ForwardSimulate <- function(f, pars, x_grid, h_grid, x0, D, z_g,
s <- x_h # also track escapement
x <- x_h # What would happen with no havest
-
+ p <- numeric(OptTime)
+
## Simulate through time ##
for(t in 1:(OptTime-1)){
@@ -43,12 +46,13 @@ ForwardSimulate <- function(f, pars, x_grid, h_grid, x0, D, z_g,
z <- z_g()
# population grows
x_h[t+1] <- z * f(x_h[t], h[t], pars) # with havest
- s[t] <- x_h[t] - q_t # anticipated escapement
+ s[t+1] <- x_h[t] - q_t # anticipated escapement
x[t+1] <- z * f(x[t], 0, pars) # havest-free dynamics
- }
+ p[t] <- profit(x_h[t], h[t])
+}
# formats output
data.frame(time = 1:OptTime, fishstock = x_h, harvest = h,
- unharvested = x, escapement = s)
+ unharvested = x, escapement = s, profit = p)
}
View
75 inst/examples/Reed.md
@@ -213,7 +213,7 @@ ggplot(subset(dt,reps==1)) +
geom_line(aes(time, harvest), col="darkgreen")
```
-![plot of chunk p0](http://farm8.staticflickr.com/7047/6983802282_7bceb5cd87_o.png)
+![plot of chunk p0](http://farm9.staticflickr.com/8015/7178088946_4ab62a1ed1_o.png)
@@ -224,10 +224,27 @@ This plot summarizes the stock dynamics by visualizing the replicates. Reed's S
```r
p1 <- ggplot(dt) + geom_abline(intercept=opt$S, slope = 0) +
geom_abline(intercept=xT, slope = 0, lty=2)
+```
+
+
+
+```
+Error: object 'xT' not found
+```
+
+
+
+```r
p1 + geom_line(aes(time, fishstock, group = reps), alpha = 0.2)
```
-![plot of chunk p1](http://farm8.staticflickr.com/7193/7129886257_51bb453c13_o.png)
+
+
+```
+Error: object 'p1' not found
+```
+
+
We can also look at the harvest dynamics:
@@ -238,7 +255,13 @@ We can also look at the harvest dynamics:
p1 + geom_line(aes(time, harvest, group = reps), alpha = 0.1, col="darkgreen")
```
-![plot of chunk p2](http://farm9.staticflickr.com/8167/6983803634_7699492b79_o.png)
+
+
+```
+Error: object 'p1' not found
+```
+
+
This strategy is supposed to be a constant-escapement strategy. We can visualize the escapement:
@@ -249,7 +272,13 @@ This strategy is supposed to be a constant-escapement strategy. We can visualize
p1 + geom_line(aes(time, escapement, group = reps), alpha = 0.1, col="darkgrey")
```
-![plot of chunk p3](http://farm8.staticflickr.com/7050/7129887427_4f690c5d65_o.png)
+
+
+```
+Error: object 'p1' not found
+```
+
+
@@ -272,10 +301,27 @@ p5 <- ggplot(policy_zoom) +
scale_colour_gradientn(colours = rainbow(4)) +
geom_abline(intercept=opt$S, slope = 0) +
geom_abline(intercept=xT, slope=0, lty=2)
+```
+
+
+
+```
+Error: object 'xT' not found
+```
+
+
+
+```r
p5
```
-![plot of chunk policy](http://farm8.staticflickr.com/7188/7129887821_5430106a9f_o.png)
+
+
+```
+Error: object 'p5' not found
+```
+
+
The harvest intensity is limited by the stock size.
@@ -290,9 +336,26 @@ p6 <- ggplot(policy_zoom) +
scale_colour_gradientn(colours = rainbow(4)) +
geom_abline(intercept=opt$S, slope = 0) +
geom_abline(intercept=xT, slope=0, lty=2)
+```
+
+
+
+```
+Error: object 'xT' not found
+```
+
+
+
+```r
p6 + geom_line(aes(time, fishstock, group = reps), alpha = 0.1, data=dt)
```
-![plot of chunk policy2](http://farm8.staticflickr.com/7108/7129888199_255a56ecb2_o.png)
+
+
+```
+Error: object 'p6' not found
+```
+
+
View
2  inst/examples/Reed_knit_.md
@@ -1,4 +1,4 @@
-`ro cache=TRUE, tidy=FALSE, warning=FALSE, comment=NA, message=FALSE, refresh=1 or`
+`ro cache=TRUE, tidy=FALSE, warning=FALSE, comment=NA, message=FALSE, refresh=2 or`
``` {r echo=FALSE }
opts_knit$set(upload.fun = socialR::flickr.url)
View
294 inst/examples/value_of_information.Rmd
@@ -0,0 +1,294 @@
+`ro cache=TRUE, tidy=FALSE, warning=FALSE, comment=NA, message=FALSE, refresh=2 fig.width=10, fig.height=10 or`
+
+``` {r echo=FALSE, cache=FALSE }
+render_jekyll()
+opts_knit$set(upload.fun = socialR::flickr.url)
+options(device = function(width = 5, height = 5) {
+ pdf(NULL, width = width, height = height)
+})
+knitcitations::cleanbib()
+````
+
+# Calculating the value of information
+ * author Carl Boettiger, <cboettig@gmail.com>
+ * license: CC0
+
+ Implements a numerical version of the SDP described in `r knitcitations::citep("10.1016/j.jeem.2004.11.005")`.
+ Compute the optimal solution under different forms of uncertainty and compare the results.
+
+``` {r setup, echo=FALSE}
+rm(list=ls())
+require(pdgControl)
+require(reshape2)
+require(ggplot2)
+require(data.table)
+rm(list=ls())
+````
+
+
+Define noise parameters
+
+
+Chose the state equation / population dynamics function
+
+``` {r }
+f <- BevHolt
+````
+
+Note that the `pdg_control` pacakge already has a definition for the `BevHolt` function, (typing the function name prints the function)
+
+``` {r }
+BevHolt
+````
+
+That is, \\( f(x,h) = \frac{A x}{1 + B x} \\)
+
+Of course we could pass in any custom function of stocksize `x`, harvest `h` and parameter vector `p` in place of `BevHolt`. Note that we would need to write this function explicitly so that it can take vector values of `x` (i.e. uses `sapply`), an annoying feature of `R` for users comming from Matlab.
+
+
+We must now define parameters for the function. Note that the positive stationary root of the model is given by \\( \frac{A-1}{B} \\), which we'll store for future reference as `K`.
+
+``` {r }
+pars <- c(1.5, 0.05)
+K <- (pars[1] - 1)/pars[2]
+````
+
+
+
+and we use a harvest-based profit function with default parameters
+
+``` {r profit}
+profit <- profit_harvest(price=1, c0 = 0.01)
+````
+
+The `profit_harvest` function has the form \\( \Pi = h - \left( c_0 + c_1 \frac{h}{x} \right) \frac{h}{x} \\), conditioned on \\( h > x \\) and \\(x > 0 \\). Note that the R code defines a function from another function using a trick known as a _closure_. Again we could write a custom profit function as long as it can take a vector stock size `x` and a scalar harvest level `h`. Details for provided functions can be found in the manual, i.e. `?profit_harvest`.
+
+
+Now we must set up the discrete grids for stock size and havest levels (which will use same resolution as for stock), in order to calculate the SDP solution. Here we set the gridsize to 100.
+
+``` {r grid}
+x_grid <- seq(0, 2 * K, length = 100)
+h_grid <- x_grid
+````
+
+
+# Scenarios:
+
+We calculate the stochastic transition matrix for the probability of going from any state \\(x_t \\) to any other state \\(x_{t+1}\\) the following year, for each possible choice of harvest \\( h_t \\). This provides a look-up table for the dynamic programming calculations.
+
+
+In the Sethi case, computing the distribution over multiple sources of noise is actually quite difficult. Simulation turns out to be more efficient than numerically integrating over each distribution.
+
+
+## No Uncertainty
+
+
+``` {r noisefns}
+sigma_g <- 0.0 # Noise in population growth
+sigma_m <- 0.0 # noise in stock assessment measurement
+sigma_i <- 0.0 # noise in implementation of the quota
+z_g <- function() rlnorm(1, 0, sigma_g) # mean 1
+z_m <- function() rlnorm(1, 0, sigma_m) # mean 1
+z_i <- function() rlnorm(1, 0, sigma_i) # mean 1
+````
+Find the transition matrix
+
+``` {r SDP_Mat}
+SDP_Mat <- determine_SDP_matrix(f, pars, x_grid, h_grid, sigma_g )
+````
+
+Find the optimum solution
+
+``` {r find_dp_opt }
+opt <- find_dp_optim(SDP_Mat, x_grid, h_grid, OptTime=25, xT=0,
+ profit, delta=0.05, reward=0)
+````
+
+Simulate
+
+``` {r simulate }
+sims_known <- lapply(1:100, function(i){
+ ForwardSimulate(f, pars, x_grid, h_grid, x0=K, opt$D, z_g, z_m, z_i, profit)
+})
+````
+
+## Growth uncertainty
+
+
+``` {r }
+sigma_g <- 0.1 # Noise in population growth
+sigma_m <- 0.0 # noise in stock assessment measurement
+sigma_i <- 0.0 # noise in implementation of the quota
+z_g <- function() rlnorm(1, 0, sigma_g) # mean 1
+z_m <- function() rlnorm(1, 0, sigma_m) # mean 1
+z_i <- function() rlnorm(1, 0, sigma_i) # mean 1
+````
+Find the transition matrix
+
+``` {r }
+SDP_Mat <- determine_SDP_matrix(f, pars, x_grid, h_grid, sigma_g )
+````
+
+Find the optimum solution
+
+``` {r }
+opt <- find_dp_optim(SDP_Mat, x_grid, h_grid, OptTime=25, xT=0,
+ profit, delta=0.05, reward=0)
+````
+
+Simulate
+
+``` {r }
+sims_g <- lapply(1:100, function(i){
+ ForwardSimulate(f, pars, x_grid, h_grid, x0=K, opt$D, z_g, z_m, z_i, profit)
+})
+````
+
+
+## Growth & stock measurement uncertainty
+
+
+``` {r }
+sigma_g <- 0.1 # Noise in population growth
+sigma_m <- 0.1 # noise in stock assessment measurement
+sigma_i <- 0.0 # noise in implementation of the quota
+z_g <- function() rlnorm(1, 0, sigma_g) # mean 1
+z_m <- function() rlnorm(1, 0, sigma_m) # mean 1
+z_i <- function() rlnorm(1, 0, sigma_i) # mean 1
+````
+
+Find the transition matrix. Use the simulation method to account for the extra uncertainties
+
+``` {r }
+require(snowfall)
+sfInit(parallel=TRUE, cpu=4)
+SDP_Mat <- SDP_by_simulation(f, pars, x_grid, h_grid, z_g, z_m, z_i, reps=999)
+````
+
+Find the optimum solution
+
+``` {r }
+opt <- find_dp_optim(SDP_Mat, x_grid, h_grid, OptTime=25, xT=0,
+ profit, delta=0.05, reward=0)
+````
+
+Simulate
+
+``` {r }
+sims_gm <- lapply(1:100, function(i){
+ ForwardSimulate(f, pars, x_grid, h_grid, x0=K, opt$D, z_g, z_m, z_i, profit)
+})
+````
+
+
+
+
+## Growth, stock measurement & implementation uncertainty
+
+
+``` {r }
+sigma_g <- 0.1 # Noise in population growth
+sigma_m <- 0.1 # noise in stock assessment measurement
+sigma_i <- 0.1 # noise in implementation of the quota
+z_g <- function() rlnorm(1, 0, sigma_g) # mean 1
+z_m <- function() rlnorm(1, 0, sigma_m) # mean 1
+z_i <- function() rlnorm(1, 0, sigma_i) # mean 1
+````
+
+Find the transition matrix. Use the simulation method to account for the extra uncertainties
+
+``` {r }
+require(snowfall)
+sfInit(parallel=TRUE, cpu=4)
+SDP_Mat <- SDP_by_simulation(f, pars, x_grid, h_grid, z_g, z_m, z_i, reps=999)
+````
+
+Find the optimum solution
+
+``` {r }
+opt <- find_dp_optim(SDP_Mat, x_grid, h_grid, OptTime=25, xT=0,
+ profit, delta=0.05, reward=0)
+````
+
+Simulate
+
+``` {r }
+sims_gmi <- lapply(1:100, function(i){
+ ForwardSimulate(f, pars, x_grid, h_grid, x0=K, opt$D, z_g, z_m, z_i, profit)
+})
+````
+
+
+
+## Summarize and plot the results
+
+R makes it easy to work with this big replicate data set. We make data tidy (melt), fast (data.tables), and nicely labeled.
+
+``` {r tidy}
+
+sims <- list(known = sims_known, growth = sims_g, growth_stock = sims_gm, growth_stock_harvest = sims_gmi)
+
+dat <- melt(sims, id=names(sims_known[[1]]))
+dt <- data.table(dat)
+setnames(dt, c("L2", "L1"), c("reps", "uncertainty")) # names are nice
+````
+
+### Plots
+
+Let's begin by looking at the dynamics of a single replicate. The line shows Reed's S, the level above which the stock should be harvested (where catch should be the difference between stock and S). To confirm that this policy is being followed, note that harvesting only occurs when the stock is above this line, and harvest is proportional to the amount by which it is above. Change the replicate `reps==` to see the results from a different replicate.
+
+``` {r onerep }
+ggplot(subset(dt,reps==1)) +
+ geom_line(aes(time, fishstock)) +
+ geom_abline(intercept=opt$S, slope = 0) +
+ geom_line(aes(time, harvest), col="darkgreen") +
+ facet_wrap(~uncertainty)
+````
+
+
+This plot summarizes the stock dynamics by visualizing the replicates. Reed's S shown again.
+
+``` {r all}
+p1 <- ggplot(dt) + geom_abline(intercept=opt$S, slope = 0)
+p1 + geom_line(aes(time, fishstock, group = reps), alpha = 0.2) + facet_wrap(~uncertainty)
+````
+
+We can also look at the harvest dynamics:
+
+``` {r harvestplot}
+p1 + geom_line(aes(time, harvest, group = reps), alpha = 0.1, col="darkgreen") + facet_wrap(~uncertainty)
+````
+
+This strategy is supposed to be a constant-escapement strategy. We can visualize the escapement:
+
+``` {r escapement}
+p1 + geom_line(aes(time, escapement, group = reps), alpha = 0.1, col="darkgrey") + facet_wrap(~uncertainty)
+````
+
+
+``` {r }
+ggplot(subset(dt,reps==1)) +
+ geom_line(aes(time, profit_fishing)) +
+ geom_line(aes(time, policy_cost), col="darkblue") + facet_wrap(~uncertainty)
+
+````
+
+``` {r }
+profits <-dt[ , sum(profit), by=c("reps", "uncertainty")]
+ggplot(profits) + geom_histogram(aes(V1)) + facet_wrap(~uncertainty)
+````
+
+Summary stats
+
+``` {r }
+profits[, mean(V1), by=uncertainty]
+profits[, sd(V1), by=uncertainty]
+````
+
+
+
+# References
+
+``` {r biblio, echo=FALSE, results="asis"}
+I(knitcitations::bibliography())
+````
View
5 man/ForwardSimulate.Rd
@@ -3,7 +3,8 @@
\title{Forward simulate given the optimal havesting policy, D}
\usage{
ForwardSimulate(f, pars, x_grid, h_grid, x0, D, z_g,
- z_m = function(x) 1, z_i = function(x) 1)
+ z_m = function(x) 1, z_i = function(x) 1,
+ profit = NULL)
}
\arguments{
\item{f}{the growth function of the escapement population
@@ -31,6 +32,8 @@
\item{z_i}{a function which returns a random number from
a distribution for the implementation uncertainty in
quotas}
+
+ \item{profit}{function, if known}
}
\value{
a data frame with the time, fishstock, harvested amount,
Please sign in to comment.
Something went wrong with that request. Please try again.