Skip to content

Commit

Permalink
prior/posterior plots added to the BH-Ricker example. Comparable but …
Browse files Browse the repository at this point in the history 
…non-optimal performance by GP.
  • Loading branch information
@cboettig
cboettig committed Dec 20, 2012
1 parent 19c32b3 commit f1711a1
Show file tree
Hide file tree
Showing 2 changed files with 116 additions and 18 deletions.
11 changes: 11 additions & 0 deletions inst/examples/reed-example.Rmd
View file
Expand Up @@ -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"}
Expand All @@ -62,6 +68,11 @@ We fit a Gaussian process with
```



```{r gp-posteriors}
```


## The transition matrix of the inferred process

```{r persistence-test}
Expand Down
123 changes: 105 additions & 18 deletions inst/examples/reed-example.md
View file
Expand Up @@ -43,15 +43,15 @@ 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))
x[t+1] = z_g(sigma_g) * f(x[t], h=0, p=p)
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
Expand All @@ -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)
```


Expand Down Expand Up @@ -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)
```



Expand All @@ -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)



Expand All @@ -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)



Expand All @@ -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)



Expand All @@ -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)



Expand Down Expand Up @@ -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)



Expand Down Expand Up @@ -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)



Expand All @@ -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)



Expand All @@ -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
```

Expand Down

0 comments on commit f1711a1

Please sign in to comment.