Admb Basic Example

After a few iterations I have a working minimal example (below). Hoping that ADMB is a bit more robust than vanilla optim out of R as I loop over data sets for the sensitivity analysis (#32). Does not seem to hold in simple example here, not sure why.

Learning ADMB

Plotting and knitr options, (can generally be ignored)

Model and parameters

f <- function(x,h,p)  x * exp(p[1] * (1 - x / p[2]) * (x - p[3]) / p[2] ) 
p <- c(1, 10, 5)
K <- 10  # approx, a li'l' less
Xo <- 6 # approx, a li'l' less

Various parameters defining noise dynamics, grid, and policy costs.

sigma_g <- 0.1
z_g <- function() rlnorm(1,0, sigma_g)
x_grid <- seq(0, 1.5 * K, length=50)
Tobs <- 40
set.seed(123)

Sample Data

x <- numeric(Tobs)
x[1] <- Xo
for(t in 1:(Tobs-1))
  x[t+1] = z_g() * f(x[t], h=0, p=p)
qplot(1:Tobs, x)
plot of chunk obs
plot of chunk obs

Maximum Likelihood “by hand”

STABLIZE = 1e-10
n = length(x)
mloglik <- function(pars){ 
  r = pars[1]; k = pars[2]; c = pars[3]; s = pars[4];
  mu = (x+STABLIZE) * exp( r * (1 - x / (k+STABLIZE)) * (x - c) / (k + STABLIZE));
  mu = pmin(1e100, mu) # avoid infinite values 
  f = 0.5 * n * log(2 * pi) + n * log(s + STABLIZE) + 0.5 * sum(x - mu + STABLIZE)^2/ (s + STABLIZE)^2;

  f
  }

Starting from the true values we mostly just shrink the noise parameter:

init <- c(p, sigma_g)
mloglik(init) #true minus loglik
[1] -35.72
o <- optim(init, mloglik, method="L", lower=1e-5, upper=1e5)
o$value
[1] -247.6
o$par
[1] 0.9967297 9.9813554 5.1742699 0.0008183

While starting from arbitrary values we still find the optim.

init <- c(1,1,1,1)  
init <- c(p, sigma_g)
mloglik(init) #true minus loglik
[1] -35.72
o <- optim(init, mloglik, method="L", lower=1e-5, upper=1e5)
o$value
[1] -247.6
o$par
[1] 0.9967297 9.9813554 5.1742699 0.0008183

Okay, now lets try admb. We use R2admb which is just a convenient way to write our data and parameters into an admb file.

# install_github("R2admb", "bbolker", subdir="R2admb") # dev version
library(R2admb)

ADMB definition

We still need to define the model using ADMB notation in the procedure section. This is mostly like R or C++, with the exception of special functions like square in place of ^2, norm2 for the sum of squares, and elem_prod istead of * for the element-wise product of two arrays. The constant pi is given as M_PI, as typical of C/C++ libraries. Where these other functions are defined I’m not sure, but some useful guides to ADMB vector/matrix operations or an (undefined) list of keywords…

The equivalent model

model <- 
paste("
PARAMETER_SECTION
  vector mu(1,n) // per capita mort prob
      
PROCEDURE_SECTION
  mu = log(x) + r * elem_prod((1 - x / k), (x - c) / k);
  f = 0.5 * n * log(2 * M_PI) + n * log(s) + 0.5 * norm2(x - exp(mu)) / square(s);
")
writeLines(model, "model.tpl")

Without explicit handling of the overflow errors, ADMB does not give us reliable estimates with arbitrary starting conditions

setup_admb("/var/admb")
[1] "/var/admb"

df <- data.frame(x=x)
params <- list(r = 1, k = 1, c = 1, s = 1) ## starting parameters
bounds <- list(r = c(1e-10, 1e3), k=c(1e-10, 1e3), c=c(1e-10, 1e3), s = c(1e-5,1e3)) ## bounds
dat <- c(list(n = nrow(df)), df)
m1 <- do_admb("model",
              data = dat,
              params = params,
              bounds = bounds,
              run.opts = run.control(checkparam="write",
                                     checkdata="write", clean=FALSE))
m1
Model file: model_gen 
Negative log-likelihood: 146.5 
Coefficients:
     r      k      c      s 
0.9992 1.0023 1.0004 9.4369

But does fine with good starting values. Hmm.. thought that was supposed to be the other way around…

params <- list(r = 1, k = 10, c = 5, s = .1) ## starting parameters

m1 <- do_admb("model",
              data = dat,
              params = params,
              bounds = bounds,
              run.opts = run.control(checkparam="write",
                                     checkdata="write"))
Warning: running command './model_gen > model_gen.out' had status 1
m1
Model file: model_gen 
Negative log-likelihood: -423.8 
Coefficients:
        r         k         c         s 
2.025e-10 9.498e+00 1.051e+01 1.000e-05

Which finds a better optim (though substantailly overfit in reality)

Hans suggests adding an error term in the function definitions rather than in limiting the bounds or log transforming the variables:

The most common plase where this goes wrong is 1/0, log(0), sqrt(0), pow(0,1) etc. Your suggestion is OK, but usually I prefer to put in log(1e-10+my_expression), sqrt(1e-10+my_expression), pow(1e-10+my_expression,1)

Misc

  • Finished off posts regarding DOIs and digital archiving. A shorter version appears in my answer to opendata.stackexchange on blog archiving.

  • Interesting discussion on using PURL redirects with SHA hash to link to repositories. For instance, the commit matching our arXiv submission of the ews-review paper is found at cboettig/ews-review/33dfb58e24ceb5861dad7cf756cffe2c5d66a655. If the hash shown in the PURL doesn’t match the one at the repository then we know something has gone wrong, since it is impossible to change the contents without changing the hash. This gives us a version-stable identifier that can always be remapped to this commit, even if Github should disappear. Of course nothing guarentees that the PURL maintainer / package author does this updating, but in principle the system is more robust than simply linking to Github.

  • In other news, I should add some automatic link checking, #94 to the notebook.