Nothing
## this returns the model to be used
## no. of parameters in the model
## do not allow to fix all params (only p-1) ??? ##Correct##
ricker <-
function(obs.error="none", fixed)
{
## Ricker model w/o obs. error
## p = 3
## model parameters: (a , b , sigma) ~ {R, R+, R+}
## parameters to monitor for convergence: (a , b , lnsigma)
cm_lik_0 <- c(
"model {",
" for (i in 1:kk) {",
" N[1,i] <- exp(x[1,i])",
" for (j in 2:T) {",
" x[j,i] ~ dnorm(mu[j,i], prcx)",
" mu[j,i] <- a + b * N[j-1,i] + x[j-1,i]",
" N[j,i] <- min(exp(x[j,i]), 10000)",
" }",
" }")
cm_end_0 <- c(
a= " a ~ dnorm(1, 0.01)",
b= " b ~ dnorm(0, 10)",
sigma= " sigma <- exp(lnsigma)",
lnsigma=" lnsigma ~ dnorm(0, 1)",
prcx= " prcx <- 1/sigma^2",
"}")
## Ricker model w/ Normal obs. error
## p = 4
## model parameters: (a , b , sigma, tau) ~ {R, R, R+, R+}
## parameters to monitor for convergence: (a , b , lnsigma, lntau)
cm_lik_N <- c(
"model {",
" for (i in 1:kk) {",
" N[1,i] <- exp(y[1,i])",
" x[1,i] <- y[1,i]",
" for (j in 2:T) {",
" x[j,i] ~ dnorm(mu[j,i], prcx)",
" mu[j,i] <- a + b * N[j-1,i] + x[j-1,i]",
" N[j,i] <- min(exp(x[j,i]), 10000)",
" y[j,i] ~ dnorm(x[j,i], prcy)",
" }",
" }")
cm_end_N <- c(
sigma= " sigma <- exp(lnsigma)",
tau= " tau <- exp(lntau)",
lnsigma=" lnsigma ~ dnorm(0, 1)",
lntau= " lntau ~ dnorm(0, 1)",
a= " a ~ dnorm(0, 0.01)",
b= " b ~ dnorm(0, 10)",
prcx= " prcx <- 1/sigma^2",
prcy= " prcy <- 1/tau^2",
"}")
## Ricker model w/ Poisson obs. error
## p = 3
## model parameters: (a , b , sigma) ~ {R, R, R+}
## parameters to monitor for convergence: (a , b , lnsigma)
cm_lik_P <- c(
"model {",
" for (i in 1:kk) {",
" N[1,i] <- O[1,i]",
" x[1,i] <- log(O[1,i])",
" for (j in 2:T) {",
" x[j,i] ~ dnorm(mu[j,i], prcx)",
" mu[j,i] <- a + b * N[j-1,i] + x[j-1,i]",
" N[j,i] <- min(exp(x[j,i]), 10000)",
" O[j,i] ~ dpois(N[j,i])",
" }",
" }")
cm_end_P <- cm_end_0
## match observation error type
obs.error <- match.arg(tolower(obs.error),
c("none", "normal", "poisson"))
## put together the model file
cm_lik <- switch(obs.error,
"none" = cm_lik_0,
"poisson" = cm_lik_P,
"normal" = cm_lik_N)
cm_end <- switch(obs.error,
"none" = cm_end_0,
"poisson" = cm_end_P,
"normal" = cm_end_N)
## number of model parameters
p <- switch(obs.error,
"none" = 3,
"poisson" = 3,
"normal" = 4)
## range of support for parameters
support <- switch(obs.error,
"none" = rbind(a=c(-Inf, Inf), b=c(-Inf, Inf),
sigma=c(.Machine$double.eps, Inf)),
"poisson" = rbind(a=c(-Inf, Inf), b=c(-Inf, Inf),
sigma=c(.Machine$double.eps, Inf)),
"normal" = rbind(a=c(-Inf, Inf), b=c(-Inf, Inf),
sigma=c(.Machine$double.eps, Inf), tau=c(.Machine$double.eps, Inf)))
colnames(support) <- c("Min", "Max")
## check range of support and put in fixed values
if (!missing(fixed)) {
if (!is.list(fixed) && is.character(fixed) && length(fixed)>1)
stop("non-numeric vector provided for fixed")
fixed <- lapply(fixed, function(z) z[1])
pp <- c("a","b","sigma")
if (all(pp %in% names(fixed)))
warning("Fixing all parameters can be a bad idea, think twice!")
if (obs.error == "normal")
pp <- c(pp, "tau")
tmp <- names(fixed)[!(names(fixed) %in% pp)]
if (length(tmp))
stop("fixed parameter ", tmp, " not found in model")
if ("a" %in% names(fixed)) {
if (is.character(fixed[["a"]])) {
cm_end["a"] <- fixed[["a"]]
} else {
if (fixed[["a"]] < support["a","Min"] || fixed[["a"]] > support["a","Max"])
stop("support for fixed parameter 'a' ill-defined")
cm_end["a"] <- paste(" a <-", round(fixed[["a"]], 4))
}
}
if ("b" %in% names(fixed)) {
if (is.character(fixed[["b"]])) {
cm_end["b"] <- fixed[["b"]]
} else {
if (fixed[["b"]] < support["b","Min"] || fixed[["b"]] > support["b","Max"])
stop("support for fixed parameter 'b' ill-defined")
cm_end["b"] <- paste(" b <-", round(fixed[["b"]], 4))
}
}
if ("sigma" %in% names(fixed)) {
if (is.character(fixed[["sigma"]])) {
cm_end["sigma"] <- fixed[["sigma"]]
} else {
if (fixed[["sigma"]] < support["sigma","Min"] || fixed[["sigma"]] > support["sigma","Max"])
stop("support for fixed parameter 'sigma' ill-defined")
cm_end["sigma"] <- paste(" sigma <-", round(fixed[["sigma"]], 4))
cm_end["lnsigma"] <- " lnsigma <- log(sigma)"
}
}
if ("tau" %in% names(fixed)) {
if (is.character(fixed[["tau"]])) {
cm_end["tau"] <- fixed[["tau"]]
} else {
if (fixed[["tau"]] < support["tau","Min"] || fixed[["tau"]] > support["tau","Max"])
stop("support for fixed parameter 'tau' ill-defined")
cm_end["tau"] <- paste(" tau <-", round(fixed[["tau"]], 4))
cm_end["lntau"] <- " lntau <- log(tau)"
}
}
fixed <- unlist(lapply(fixed, function(z)
ifelse(is.numeric(z), z, NA)))
fixed <- fixed[!is.na(fixed)]
if (!length(fixed))
fixed <- NULL
} else fixed <- NULL
## put together stuff
model <- structure(c(cm_lik, unname(cm_end)),
class = "custommodel")
fancy <- c("Ricker",
ifelse(obs.error=="none", NA,
gsub("\\b(\\w)", "\\U\\1", obs.error, perl=TRUE)))
## list params (fixed will not be a parameter)
params <- switch(obs.error,
"none" = c("a","b","lnsigma"),
"poisson" = c("a","b","lnsigma"),
"normal" = c("a","b","lnsigma","lntau"))
varnames <- switch(obs.error,
"none" = c("a","b","sigma"),
"poisson" = c("a","b","sigma"),
"normal" = c("a","b","sigma","tau"))
params <- params[!(varnames %in% names(fixed))]
## this scales diagnostic parameters to the scale of the summaries
backtransf <- function(mcmc, obs.error) {
#mcmc <- as.mcmc.list(mcmc)
mcmc <- as(mcmc, "MCMClist")
vn <- varnames(mcmc)
for (i in seq_len(nchain(mcmc))) {
if ("lnsigma" %in% vn)
mcmc[[i]][,"lnsigma"] <- exp(mcmc[[i]][,"lnsigma"])
if ("lntau" %in% vn)
mcmc[[i]][,"lntau"] <- exp(mcmc[[i]][,"lntau"])
}
if ("lnsigma" %in% vn)
vn[vn=="lnsigma"] <- "sigma"
if ("lntau" %in% vn)
vn[vn=="lntau"] <- "tau"
varnames(mcmc) <- vn
mcmc
}
## this scales summaries to the scale of diagnostic parameters
transf <- function(mcmc, obs.error) {
#mcmc <- as.mcmc.list(mcmc)
mcmc <- as(mcmc, "MCMClist")
vn <- varnames(mcmc)
for (i in seq_len(nchain(mcmc))) {
if ("sigma" %in% vn)
mcmc[[i]][,"sigma"] <- log(mcmc[[i]][,"sigma"])
if ("tau" %in% vn)
mcmc[[i]][,"tau"] <- log(mcmc[[i]][,"tau"])
}
if ("sigma" %in% vn)
vn[vn=="sigma"] <- "lnsigma"
if ("tau" %in% vn)
vn[vn=="tau"] <- "lntau"
varnames(mcmc) <- vn
mcmc
}
dfun <- switch(obs.error,
## data on the log scale (logx)
## w/o missing data, log is the log data
## w/ missing values both logx (missing) and data should be provided
"none" = function(logx, mle, data,
null_obserror=FALSE, alt_obserror=FALSE) {
T <- length(logx)
m <- which(is.na(data))
if (length(m) > 0) {
if (alt_obserror) {
## null is NOE, alt is OE, NA present (II.a)
stop("not yet implemented")
} else {
## null is NOE, alt is NOE, NA present (II.b)
y <- data
y[m] <- logx[m]
rval <- sum(dnorm(y[-1],
mean = y[-T] + mle["a"] + mle["b"] * exp(y[-T]),
sd = mle["sigma"], log=TRUE))
}
} else {
logd1 <- dnorm(data[-1],
mean = data[-T] + mle["a"] + mle["b"] * exp(data[-T]),
sd = mle["sigma"], log=TRUE)
## null is NOE, alt is OE, NA absent (III.a)
logd2 <- if (alt_obserror) {
dnorm(logx[-1],
mean = logx[-T] + mle["a"] + mle["b"] * exp(logx[-T]),
sd = mle["sigma"], log=TRUE)
## null is NOE, alt is NOE, NA absent (III.b)
} else {
# dnorm((logx[-1])[m-1],
# mean = (logx[-T])[m] + mle["a"] +
# mle["b"] * exp((logx[-T])[m]),
# sd = mle["sigma"], log=TRUE)
0
}
rval <- sum(logd1) + sum(logd2)
}
rval
},
## data on the log scale
"poisson" = function(logx, mle, data,
null_obserror=FALSE, alt_obserror=FALSE) {
T <- length(data)
if (!(!null_obserror && any(is.na(data)))) {
logd1 <- dnorm(logx[-1],
mean = logx[-T] + mle["a"] + mle["b"] * exp(logx[-T]),
sd = mle["sigma"], log=TRUE)
logd2 <- dpois(data[-1], exp(logx[-1]), log=TRUE)
rval <- sum(logd1) + sum(logd2, na.rm=TRUE)
} else {
stop("not yet implemented")
}
rval
},
## data on the log scale
"normal" = function(logx, mle, data,
null_obserror=FALSE, alt_obserror=FALSE) {
T <- length(data)
if (!(!null_obserror && any(is.na(data)))) {
logd1 <- dnorm(logx[-1],
mean = logx[-T] + mle["a"] + mle["b"] * exp(logx[-T]),
sd = mle["sigma"], log=TRUE)
logd2 <- dnorm(data[-1],
mean = logx[-1],
sd = mle["tau"], log=TRUE)
rval <- sum(logd1) + sum(logd2, na.rm=TRUE)
} else {
stop("not yet implemented")
}
rval
})
neff <- function(obs)
sum(!is.na(obs)) - 1
out <- new("pvamodel")
out@growth.model <- "ricker"
out@obs.error <- obs.error
out@model <- model
out@p <- as.integer(p)
out@support <- support
out@params <- params
out@varnames <- varnames
out@fixed <- fixed
out@fancy <- fancy
out@backtransf <- backtransf
out@transf <- transf
out@logdensity <- dfun
out@neffective <- neff
out
}
## make this as an S4 class??? use contains???
#ricker()
#ricker("Pois")
#ricker("normal")
#ricker("normal", fixed=c(a=5, sigma=0.5))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.