diff --git a/inst/examples/reed-example.Rmd b/inst/examples/reed-example.Rmd index 54d8070..d09b864 100644 --- a/inst/examples/reed-example.Rmd +++ b/inst/examples/reed-example.Rmd @@ -48,6 +48,12 @@ Simulate data ```{r par-est} ``` + + +```{r gp-priors} +``` + + Estimates a Ricker curve with parameters $r =$ `r p_alt[1]` and $K =$ `r p_alt[2]` ```{r gp-fit, dependson="lag-data"} @@ -62,6 +68,11 @@ We fit a Gaussian process with ``` + +```{r gp-posteriors} +``` + + ## The transition matrix of the inferred process ```{r persistence-test} diff --git a/inst/examples/reed-example.md b/inst/examples/reed-example.md index 462cc75..947cd09 100644 --- a/inst/examples/reed-example.md +++ b/inst/examples/reed-example.md @@ -43,7 +43,7 @@ We consider stochastic growth driven by a lognormal noise process, $X_{t+1} = z_ ```r -Tobs <- 30 +Tobs <- 50 x <- numeric(Tobs) x[1] <- x_0_observed for(t in 1:(Tobs-1)) @@ -51,7 +51,7 @@ for(t in 1:(Tobs-1)) plot(x) ``` -![plot of chunk sim-obs](http://carlboettiger.info/assets/figures/2012-12-18-14275a382d-sim-obs.png) +![plot of chunk sim-obs](http://carlboettiger.info/assets/figures/2012-12-19-19c32b33a1-sim-obs.png) Simulate data @@ -75,16 +75,39 @@ sigma_g_alt <- o$par['s'] ``` -Estimates a Ricker curve with parameters $r =$ `0.3906` and $K =$ `6.3984` + + + +```r +s2.p <- c(50,50) +tau2.p <- c(20,1) +d.p = c(10, 1/0.01, 10, 1/0.01) +nug.p = c(10, 1/0.01, 10, 1/0.01) +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) +``` + + + +Estimates a Ricker curve with parameters $r =$ `0.4182` and $K =$ `7.7853` ```r gp <- bgp(X=obs$x, XX=x_grid, Z=obs$y, verb=0, - meanfn="linear", bprior="b0", BTE=c(2000,6000,2), m0r1=FALSE, - corr="exp", trace=TRUE, beta = c(0,0), - s2.p = c(50,50), d.p = c(10, 1/0.01, 10, 1/0.01), nug.p = c(10, 1/0.01, 10, 1/0.01), - s2.lam = "fixed", d.lam = "fixed", nug.lam = "fixed", - tau2.lam = "fixed", tau2.p = c(50,1)) + 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") +``` + +``` +Warning: for memory/storage reasons, trace not recommended when +3*(10+d)*(BTE[2]-BTE[1])*R*(nn+1)/BTE[3]=23562000 > 1e+7. Try reducing +nrow(XX) ``` @@ -116,7 +139,71 @@ ggplot(tgp_dat) + geom_ribbon(aes(x,y,ymin=ymin,ymax=ymax), fill="gray80") + scale_colour_manual(values=cbPalette) ``` -![plot of chunk gp-plot](http://carlboettiger.info/assets/figures/2012-12-18-14275a382d-gp-plot.png) +![plot of chunk gp-plot](http://carlboettiger.info/assets/figures/2012-12-19-19c32b33a1-gp-plot.png) + + + + + +```r +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_density(aes(value), lwd=2) + + 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") +``` + +``` +stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust +this. +``` + +``` +stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust +this. +``` + +``` +stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust +this. +``` + +``` +stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust +this. +``` + +``` +stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust +this. +``` + +``` +stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust +this. +``` + +![plot of chunk gp-posteriors](http://carlboettiger.info/assets/figures/2012-12-19-19c32b33a1-gp-posteriors1.png) + +```r +ggplot(prior_curves) + + geom_line(aes(x=value, y=density), col="red", lwd=2) + + facet_wrap(~ variable, scale="free") +``` + +![plot of chunk gp-posteriors](http://carlboettiger.info/assets/figures/2012-12-19-19c32b33a1-gp-posteriors2.png) + +```r +#ggplot(subset(posteriors, variable=="nug")) + geom_histogram(aes(x=value, y = ..density..), lwd=2) + stat_function(fun = nug_prior, col="red", lwd=2) +#ggplot(subset(posteriors, variable=="s2")) + geom_histogram(aes(x=value, y = ..density..), lwd=2) + stat_function(fun = s2_prior, col="red", lwd=2) +``` @@ -140,7 +227,7 @@ for(s in 1:OptTime) qplot(x_grid, xt10[1,]) + geom_point(aes(y=xt1[1,]), col="grey") ``` -![plot of chunk gp-F-sim](http://carlboettiger.info/assets/figures/2012-12-18-14275a382d-gp-F-sim.png) +![plot of chunk gp-F-sim](http://carlboettiger.info/assets/figures/2012-12-19-19c32b33a1-gp-F-sim.png) @@ -153,7 +240,7 @@ for(s in 1:OptTime) qplot(x_grid, yt10[1,]) + geom_point(aes(y=yt1[1,]), col="grey") ``` -![plot of chunk par-F-sim](http://carlboettiger.info/assets/figures/2012-12-18-14275a382d-par-F-sim.png) +![plot of chunk par-F-sim](http://carlboettiger.info/assets/figures/2012-12-19-19c32b33a1-par-F-sim.png) @@ -162,7 +249,7 @@ transition <- melt(data.frame(x = x_grid, gp = xt1[1,], parametric = yt1[1,]), i ggplot(transition) + geom_point(aes(x,value, col=variable)) ``` -![plot of chunk F-sim-plot](http://carlboettiger.info/assets/figures/2012-12-18-14275a382d-F-sim-plot.png) +![plot of chunk F-sim-plot](http://carlboettiger.info/assets/figures/2012-12-19-19c32b33a1-F-sim-plot.png) @@ -175,7 +262,7 @@ for(s in 1:OptTime) qplot(x_grid, zt10[1,]) + geom_point(aes(y=zt1[1,]), col="grey") ``` -![plot of chunk est-F-sim](http://carlboettiger.info/assets/figures/2012-12-18-14275a382d-est-F-sim.png) +![plot of chunk est-F-sim](http://carlboettiger.info/assets/figures/2012-12-19-19c32b33a1-est-F-sim.png) @@ -216,7 +303,7 @@ policy_plot <- ggplot(policies, aes(stock, stock - value, color=method)) + policy_plot ``` -![plot of chunk policy_plot](http://carlboettiger.info/assets/figures/2012-12-18-14275a382d-policy_plot.png) +![plot of chunk policy_plot](http://carlboettiger.info/assets/figures/2012-12-19-19c32b33a1-policy_plot.png) @@ -269,7 +356,7 @@ ggplot(dt) + scale_colour_manual(values=cbPalette, guide = guide_legend(override.aes = list(alpha = 1))) ``` -![plot of chunk sim-fish](http://carlboettiger.info/assets/figures/2012-12-18-14275a382d-sim-fish.png) +![plot of chunk sim-fish](http://carlboettiger.info/assets/figures/2012-12-19-19c32b33a1-sim-fish.png) @@ -280,7 +367,7 @@ ggplot(dt) + scale_colour_manual(values=cbPalette, guide = guide_legend(override.aes = list(alpha = 1))) ``` -![plot of chunk sim-harvest](http://carlboettiger.info/assets/figures/2012-12-18-14275a382d-sim-harvest.png) +![plot of chunk sim-harvest](http://carlboettiger.info/assets/figures/2012-12-19-19c32b33a1-sim-harvest.png) @@ -293,8 +380,8 @@ cbind(means, sd = sds$V1) ``` method V1 sd -1: GP 23.99 0.3890 -2: Parametric 23.59 0.4076 +1: GP 23.83 0.3527 +2: Parametric 23.65 0.4774 3: True 23.68 0.4513 ```