Nothing
#####################################################
# Performs simulated annealing to find maximum likelihood
# estimates for a set of parameters.
#
# Arguments:
# model = model function to parameterize. Arguments to this function will be
# provided from par, var, and source_data.
#
# par = list of parameters for which we are using simulated annealing to
# find maximum likelihood estimates. Each par element name matches an
# argument in a function (either model or pdf). All values listed in par
# must be numeric vectors. Vectors of length greater than one have each of
# their elements treated separately as individual parameters to estimate.
#
# var = list of other variables and data needed by the model and pdf
# functions, any type as needed. These will be kept constant.
#
# source_data = data frame with dependent variable and associated
# independent variables.
#
# par_lo = list of lower bounds for each parameter. Each element must match
# an element in par in both name and size. Ommitted values are assumed to be
# negative infinity.
#
# par_hi = list of upper bounds for each parameter. Each element must match
# an element in par in both name and size. Ommitted values are assumed to
# be positive infinity.
#
# pdf = probability density function. Make sure to use a function such as
# dnorm that can calculate log probability.
#
# dep_var = dependent variable label in source_data.
#
# initial_temp = temperature at which to start annealing.
#
# temp_red = interval by which to reduce temperature.
#
# max_iter = maximum iterations, where one iteration is one varying of each
# parameter to parameterize. The annealing will run for this many
# iterations, unless alternate conditions have been specified (see min_change
# and min_drops).
#
# dep_var = label of model dependent variable
#
# ns = interval between changes in range
#
# nt = interval between drops in temperature
#
# min_change = an alternate way to specify quitting conditions. This is the
# minimum amount of change in likelihood in min_drop number of temperature
# drops. If the change is less, execution stops.
#
# min_drops = the companion to min_change for alternate quitting conditions.
# This is the number of temperature drops over which the likelihood must
# have changed more than min_change for execution to continue.
#
# hessian = if TRUE, include the standard errors and the Hessian matrix
# in the output. If FALSE, do not.
#
# delta = when calculating support limits, the number of pieces into which
# to divide each parameter. This is the size of the "step" the function
# takes in trying to find the support limits.
#
# slimit = when calculating support limits, the number of likelihood
# units less than the optimum likelihood for which the support intervals
# will be calculated. 2 units is semi-standard. 1.92 units corresponds
# roughly to a 95% confidence interval.
#
# c = range reduction parameter
#
# note = any note to self the user wants to make. This will be
# included in the output.
#
# ... = additional arguments to model, pdf, etc. This is not a recommended
# way to pass additional arguments but I have chosen to support it.
#
# Author: Lora Murphy, Cary Institute of Ecosystem Studies
# murphyl@caryinstitute.org
######################################################
anneal<-function(model, par, var, source_data, par_lo = NULL,
par_hi = NULL, pdf, dep_var, initial_temp = 3, temp_red = 0.95,
ns = 20, nt = 100, max_iter = 50000, min_change = 0, min_drops = 100,
hessian = TRUE, delta = 100, slimit = 2, c = 2, note = "", show_display = TRUE, ...) {
##
## Error checking
##
# Check types
if (!is.function(model)) {
stop("anneal: model is not a function.\n")
}
if (!is.list(par)) {
stop("anneal: par is not a list.\n")
}
if (!is.list(var)) {
stop("anneal: var is not a list.\n")
}
if (!is.null(par_lo) && !is.list(par_lo)) {
stop("anneal: par_lo is not a list.\n")
}
if (!is.null(par_hi) && !is.list(par_hi)) {
stop("anneal: par_hi is not a list.\n")
}
if (!is.data.frame(source_data)) {
stop("anneal: source_data is not a data frame.\n")
}
if (length(source_data) < 1) {
stop("anneal: source_data contains no data.\n")
}
if (length(source_data[[1]]) < 1) {
stop("anneal: source_data contains no data.\n")
}
if (!is.function(pdf)) {
stop("anneal: pdf is not a function.\n")
}
# Make sure there is something in par
parnames <- names(par)
numpars <- length(par)
if (numpars < 1) {
stop("anneal: no parameters in par - nothing to anneal!\n")
}
# Make sure all values in par are numeric
for (i in 1:numpars) {
if (!is.numeric(par[[i]])) {
stop("anneal: All values in par must be numeric.\n")
}
if (!is.vector(par[[i]])) {
stop("anneal: All values in par must be vectors.\n")
}
}
# Make sure the dependent variable is part of source_data
if (!is.character(dep_var)) {
stop("anneal: I cannot understand the dependent variable name. I expect a string.\n")
}
if(!any(names(source_data)==dep_var)) {
stop("anneal: the dependent variable must be part of source_data.\n")
}
# Make sure that par_lo and par_hi have corresponding labels in par,
# and that the vectors are the same length
if (is.null(par_lo)) {
par_lo <- par
for (i in 1:numpars) {
par_lo[[parnames[i]]] = rep(-.Machine$double.xmax, times = length(par[[i]]))
}
}
if (is.null(par_hi)) {
par_hi <- par
for (i in 1:numpars) {
par_hi[[parnames[i]]] = rep(.Machine$double.xmax, times = length(par[[i]]))
}
}
for (i in 1:length(names(par_lo))) {
if(!any(parnames == names(par_lo)[i])) {
stop("anneal: all values in par_lo must be present in par. Can't find value ", names(par_lo)[i], ".\n")
}
if (length(par_lo[[i]]) != length( par[(parnames==names(par_lo)[i])][[1]] )) {
stop("anneal: par_lo value ", names(par_lo)[i], " must be of the same length as its corresponding element in par.\n")
}
}
for (i in 1:length(names(par_hi))) {
if(!any(parnames == names(par_hi)[i])) {
stop("anneal: all values in par_hi must be present in par. Can't find value ", names(par_hi)[i], ".\n")
}
if (length(par_hi[[i]]) != length( par[(parnames==names(par_hi)[i])][[1]] )) {
stop("anneal: par_hi value ", names(par_hi)[i], " must be of the same length as its corresponding element in par.\n")
}
}
# For each value in par, if the value in par_lo or par_hi is
# missing, replace it with "infinity"
# I used to add to the old array by name but this created arrays in different
# orders - very bad for later!
new_par_lo <- list()
new_par_hi <- list()
for (i in 1:numpars) {
if (!any(names(par_lo) == parnames[i])) {
new_par_lo[[i]] = rep(-.Machine$double.xmax, times = length(par[[i]]))
} else new_par_lo[[i]] = par_lo[[parnames[i]]]
if (!any(names(par_hi) == parnames[i])) {
new_par_hi[[i]] = rep(.Machine$double.xmax, times = length(par[[i]]))
} else new_par_hi[[i]] = par_hi[[parnames[i]]]
}
names(new_par_lo) <- parnames
names(new_par_hi) <- parnames
par_lo <- new_par_lo
par_hi <- new_par_hi
new_par_lo <- NULL
new_par_hi <- NULL
# Replace any true values of infinity with our fake infinity
for (i in 1:length(par_lo)) {
par_lo[[i]] = ifelse(is.infinite(par_lo[[i]]), -.Machine$double.xmax, par_lo[[i]])
}
for (i in 1:length(par_hi)) {
par_hi[[i]] = ifelse(is.infinite(par_hi[[i]]), .Machine$double.xmax, par_hi[[i]])
}
##
## Prep the likelihood calculations
##
# Here's where we'll put results of all pre-evaluations
eval_results <- NULL
# This is where to put the predicted
predicted <- NULL
# Combine par and var into one - we will put the varying values here
par <- c(par, var)
# Set up the PDF call - set it up to use a copy of par
datasets<-NULL
datasets[[1]]<-list(value=source_data, varname="source_data")
pdf_call<-analyze_function(pdf, list(value=par, varname="par_copy"), NULL, datasets, ...)
# Set up the model call - set it up to use a copy of par
model_call<-analyze_function(model, list(value=par, varname="par_copy"), NULL, datasets, ...)
# Initialize support limits to zeroes before we "flatten" arrays
upper_limit<-list()
lower_limit<-list()
for (i in 1:length(par_lo)) {
upper_limit[[names(par_lo)[[i]]]]<-rep(0, length(par_lo[[i]]))
lower_limit[[names(par_lo)[[i]]]]<-rep(0, length(par_lo[[i]]))
}
# If there are any vectors of values in a varying parameter, "flatten"
# into 1-d vectors; create an array of index vectors for
# pinpointing just one parameter value in par, even if it's in
# an array
#
# Plus make a flat array of the initial values for calculating step
# Save copies for output
par_lo_copy <- par_lo
par_hi_copy <- par_hi
# Make a placeholder for step
par_step_copy <- par_lo
lo_copy<-NULL
hi_copy<-NULL
initial_vals <- NULL
nm <- NULL # new names
par_index<-NULL
for (i in 1:length(par_hi)) {
numpars <- numpars + (length(par_hi[[i]]) - 1)
# Create new names for display
if (length(par_hi[[i]]) == 1) nm<-c(nm, parnames[[i]])
else nm <- c(nm, paste(parnames[[i]], c(1:length(par_hi[[i]])), sep=""))
lo_copy<-c(lo_copy,par_lo[[i]])
hi_copy<-c(hi_copy,par_hi[[i]])
initial_vals<-c(initial_vals, par[[parnames[[i]]]])
for (j in 1:length(names(par)))
if (parnames[[i]] == names(par)[[j]]) ind <- j
for (j in 1:length(par_hi[[i]])) {
par_index[[length(par_index)+1]] <- c(ind, j)
}
}
par_lo<-lo_copy
par_hi<-hi_copy
lo_copy<-NULL
hi_copy<-NULL
names(par_lo) <- nm
names(par_hi) <- nm
names(initial_vals) <- nm
nm <- NULL
##
## Since I can't break the habits of C++ - a rundown of my variables...
##
iterate <- TRUE # flag for whether or not to keep cycling
using <- 1 # index for the parameter we're working with
cycles <- 0 # how many cycles of annealing we've done
best_lh <- NULL # best likelihood
acc_lh <- NULL # last accepted likelihood
best_par <- par # set of parameters providing the best likelihood
best_cycle <- 1 # number when the current best likelihood was calculated
best_predicted<-NULL # Predicted values that go with the best parameters
nacc <- 0 # number of accepted likelihood values
temp <- initial_temp # annealing temperature
lhist <- NULL # array of lists with likelihood history
#par_copy <- NULL # copy of parameters that we will vary
#parnames <- names(par_step) # names of parameters to vary - already defined above
#numpars # number of parameters to vary - already defined above
#par_index # index array for retrieving values from par - already defined above
#upper_limit # upper support limit for all varying parameters - already defined above
#lower_limit # lower support limit for all varying parameters - already defined above
#par_step # search interval - defined below
nacp <- rep(0, numpars) # number of times we accepted changes to each parameter
range_cycles <- numpars * ns # number of loops before we adjust range
temp_cycles <- numpars * ns * nt # number of loops before we reduce temperature
display_cycles<-numpars*100 # regular update point
numobs <- length(source_data[[dep_var]]) # number of observations
aiccorr <- NULL # value for AIC corr (not entirely sure what that is)
sst <- NULL # used for calculating R2 but not very sure what it is
sse <- NULL # used for calculating R2 but not very sure what it is
R2 <- NULL # R2 (squared) value (NOT r2) - goodness-of-fit, or proportion
# of variance explained by the model
slp <- NULL # slope of observed regressed on predicted
sumx2 <- NULL # used for calculating slp
sumxy <- NULL # used for calculating slp
temp_hist_inds <- 1 # indexes in likelihood history of each temperature reduction
# Calculate the initial search interval, which for each parameter is the
# greater distance between the initial value and the upper and lower bound
par_step = numeric()
for (i in 1:length(par_hi))
par_step[i] <- max(abs(par_hi[i] - initial_vals[i]), abs(initial_vals[i] - par_lo[i]))
names(par_step) <- names(par_hi)
initial_vals<-NULL
# Calculate SST for observed values
sst <- mean(source_data[[dep_var]])
sst <- (source_data[[dep_var]] - sst)
sst <- sst * sst
sst <- sum(sst)
##
## Do the first likelihood calculation with initial values
##
par_copy <- par
# First: any sub-functions needed by the model
if (!is.null(model_call$pre_eval)) {
level<-1
pre_evals<-NULL
pre_evals[[level]]<-model_call$pre_eval
curr_index<-1
max_index<-NULL
max_index[[level]]<-length(model_call$pre_eval)
while (level > 0) {
if (is.null(pre_evals[[level]][[curr_index[[level]]]]$pre_eval)) {
eval_results[[pre_evals[[level]][[curr_index[[level]]]]$parname]]<-do.call(pre_evals[[level]][[curr_index[[level]]]]$fun, pre_evals[[level]][[curr_index[[level]]]]$call)
curr_index[[level]]<-curr_index[[level]]+1
if (curr_index[[level]] > max_index[[level]]) {
# We've evaluated all the pre-evals in this level - back up
# one level
level <- level - 1
}
}
else {
# The current level has pre-evals - before evaluating the current level,
# go down a level to the first of its pre-evals
level <- level + 1
curr_index[[level]] <- 1
max_index[[level]] <- length(pre_evals[[level]][[curr_index[[level]]]]$pre_eval)
pre_evals[[level]] <- pre_evals[[level]][[curr_index[[level]]]]$pre_eval
}
}
}
eval_results <- as.list(eval_results)
# Second: the model predicted values
predicted <- do.call(model,model_call$call)
if (any(is.infinite(predicted)) || any(is.nan(predicted)) || any(is.na(predicted))) {
stop("anneal: initial conditions caused the model to produce math errors.\n")
}
# Third: any sub-functions needed by the pdf
if (!is.null(pdf_call$pre_eval)) {
level<-1
pre_evals<-NULL
pre_evals[[level]]<-pdf_call$pre_eval
curr_index<-1
max_index<-NULL
max_index[[level]]<-length(pdf_call$pre_eval)
while (level > 0) {
if (is.null(pre_evals[[level]][[curr_index[[level]]]]$pre_eval)) {
eval_results[[pre_evals[[level]][[curr_index[[level]]]]$parname]]<-do.call(pre_evals[[level]][[curr_index[[level]]]]$fun, pre_evals[[level]][[curr_index[[level]]]]$call)
curr_index[[level]]<-curr_index[[level]]+1
if (curr_index[[level]] > max_index[[level]]) {
# We've evaluated all the pre-evals in this level - back up
# one level
level <- level - 1
}
}
else {
# The current level has pre-evals - before evaluating the current level,
# go down a level to the first of its pre-evals
level <- level + 1
curr_index[[level]] <- 1
max_index[[level]] <- length(pre_evals[[level]][[curr_index[[level]]]]$pre_eval)
pre_evals[[level]] <- pre_evals[[level]][[curr_index[[level]]]]$pre_eval
}
}
}
# Fourth: sum over the PDF
best_lh <- sum (do.call(pdf,pdf_call$call)) # best likelihood
best_predicted <- predicted
if (is.infinite(best_lh) || is.nan(best_lh) || is.na(best_lh)) best_lh <- -1000000
acc_lh <- best_lh
lhistcycles <- cycles
lhisthood <- best_lh
lhisttemp <- temp
lhistpar <- vector(mode="list")
temppar <- list()
for (q in 1:length(parnames)) temppar[[parnames[q]]] <- par[[parnames[q]]]
lhistpar[[1]] <- temppar
# Calculate slope of observed regressed on predicted
sumx2 <- sum(predicted * predicted)
sumxy <- sum(predicted * source_data[[dep_var]])
slp <- sumxy/sumx2
# Calculate R2 = 1 - (SSE/SST)
sse <- (source_data[[dep_var]] - predicted)
sse <- sse * sse
sse <- sum(sse)
R2 <- 1 - (sse/sst)
# Calculate AIC corr (corrected for small sample size)
aiccorr <- (-2.0 * best_lh)+((2*numpars)*(numobs/(numobs - numpars - 1)))
# Setup output display and show it for the first time
# Create a window so likdisplay will default to showing the graphics
if (show_display) {
plot(x=1,y=1)
layout(matrix(c(1,2), 2, 1, byrow = TRUE)) # this makes a two-part graph
likdisplay(lhistcycles, lhisthood, slp, R2, aiccorr, temp, max_iter)
}
tryCatch ({
##
## DO ANNEALING
##
# Keep cycling until we're done
while (iterate) {
cycles <- cycles + 1;
# Make a copy of parameters before we vary
par_copy <- par
# Choose a new value for the parameter being varied this time
par_copy[[par_index[[using]]]] <- par[[par_index[[using]]]] + ((runif(1)*2 - 1)* par_step[[using]]);
# Make sure the new value is in bounds
if (par_copy[[par_index[[using]]]] < par_lo[[using]]) {
par_copy[[par_index[[using]]]] <- par[[par_index[[using]]]] - (runif(1) * (par[[par_index[[using]]]]- par_lo[[using]]))
}
if (par_copy[[par_index[[using]]]] > par_hi[[using]]) {
par_copy[[par_index[[using]]]] <- par[[par_index[[using]]]] + (runif(1) * (par_hi[[using]]-par[[par_index[[using]]]]))
}
##
## Calculate likelihood with this set of parameters
##
# First: any sub-functions needed by the model
if (!is.null(model_call$pre_eval)) {
level<-1
pre_evals<-NULL
pre_evals[[level]]<-model_call$pre_eval
curr_index<-1
max_index<-NULL
max_index[[level]]<-length(model_call$pre_eval)
while (level > 0) {
if (is.null(pre_evals[[level]][[curr_index[[level]]]]$pre_eval)) {
# Some of these may actually have stored a copy of "par_copy" - refresh if this is the case
if (any(names(pre_evals[[level]][[curr_index[[level]]]]$call)=="par_copy")) {
pre_evals[[level]][[curr_index[[level]]]]$call$par_copy <- par_copy
}
eval_results[[pre_evals[[level]][[curr_index[[level]]]]$parname]]<-do.call(pre_evals[[level]][[curr_index[[level]]]]$fun, pre_evals[[level]][[curr_index[[level]]]]$call)
curr_index[[level]]<-curr_index[[level]]+1
if (curr_index[[level]] > max_index[[level]]) {
# We've evaluated all the pre-evals in this level - back up
# one level
level <- level - 1
}
}
else {
# The current level has pre-evals - before evaluating the current level,
# go down a level to the first of its pre-evals
level <- level + 1
curr_index[[level]] <- 1
max_index[[level]] <- length(pre_evals[[level]][[curr_index[[level]]]]$pre_eval)
pre_evals[[level]] <- pre_evals[[level]][[curr_index[[level]]]]$pre_eval
}
}
}
eval_results <- as.list(eval_results)
# Second: the model predicted values
predicted <- do.call(model,model_call$call)
# Make sure that all values are defined
if (!any(is.infinite(predicted)) && !any(is.nan(predicted)) && !any(is.na(predicted))) {
# Third: any sub-functions needed by the pdf
if (!is.null(pdf_call$pre_eval)) {
level<-1
pre_evals<-NULL
pre_evals[[level]]<-pdf_call$pre_eval
curr_index<-1
max_index<-NULL
max_index[[level]]<-length(pdf_call$pre_eval)
while (level > 0) {
if (is.null(pre_evals[[level]][[curr_index[[level]]]]$pre_eval)) {
eval_results[[pre_evals[[level]][[curr_index[[level]]]]$parname]]<-do.call(pre_evals[[level]][[curr_index[[level]]]]$fun, pre_evals[[level]][[curr_index[[level]]]]$call)
curr_index[[level]]<-curr_index[[level]]+1
if (curr_index[[level]] > max_index[[level]]) {
# We've evaluated all the pre-evals in this level - back up
# one level
level <- level - 1
}
}
else {
# The current level has pre-evals - before evaluating the current level,
# go down a level to the first of its pre-evals
level <- level + 1
curr_index[[level]] <- 1
max_index[[level]] <- length(pre_evals[[level]][[curr_index[[level]]]]$pre_eval)
pre_evals[[level]] <- pre_evals[[level]][[curr_index[[level]]]]$pre_eval
}
}
}
# Fourth: sum over the PDF
lh <- sum (do.call(pdf,pdf_call$call))
if (is.infinite(lh) || is.nan(lh) || is.na(lh)) lh <- -1000000
# Accept the new values if likelihood increases or at least stays the same
# since the last accepted likelihood
if (lh >= acc_lh) {
par <- par_copy
acc_lh <- lh
nacc <- nacc + 1
nacp[[using]] <- nacp[[using]] + 1
# If this is a new maximum, then update the maximum likelihood
if (acc_lh > best_lh) {
best_par <- par_copy
best_lh <- acc_lh
best_cycle <- cycles
best_predicted <- predicted
# Calculate slope of observed regressed on predicted
sumx2 <- sum(predicted * predicted)
sumxy <- sum(predicted * source_data[[dep_var]])
slp <- sumxy/sumx2
# Calculate R2 = 1 - (SSE/SST)
sse <- (source_data[[dep_var]] - predicted)
sse <- sse * sse
sse <- sum(sse)
R2 <- 1 - (sse/sst)
# Calculate AIC corrected for small sample size
aiccorr <- (-2.0 * best_lh)+((2*numpars)*(numobs/(numobs - numpars - 1)))
lhisttemp[length(lhisttemp)+1] <- temp
lhistcycles[length(lhistcycles) + 1] <- cycles / numpars
lhisthood[length(lhisthood) + 1] <- best_lh
temppar <- list()
for (q in 1:length(parnames)) temppar[[parnames[q]]] <- best_par[[parnames[q]]]
lhistpar[[length(lhistpar) + 1]] <- temppar
# Display
if (show_display) {
likdisplay(lhistcycles, lhisthood, slp, R2, aiccorr, temp, max_iter)
}
}
}
else {
# Use Metropolis criteria to determine whether to accept a downhill move
# lh < f, so the code below is a shortcut for exp(-1.0(abs(acc_lh-lh)/t)
prob <- exp((lh-acc_lh)/temp); # temp = current temperature
if (runif(1) < prob) {
par <- par_copy
acc_lh <- lh
nacp[[using]] <- nacp[[using]] + 1
}
}
# After the user-specified interval, adjust step so that half of
# evaluations are accepted
if (cycles %% range_cycles == 0) {
for (i in 1:numpars) {
if (par_step[[i]] != 0) {
ratio <- nacp[[i]] / ns
# C controls the adjustment of range - references suggest
# setting at 2.0. I've hard-coded the 2 since Charlie doesn't
# seem to allow setting it anywhere.
if (ratio > 0.6) par_step[[i]] <- par_step[[i]] * (1.0 + c*((ratio - 0.6)/0.4))
else {if (ratio < 0.4) par_step[[i]] <- par_step[[i]] / (1.0+c*((0.4 - ratio)/0.4))}
if (is.infinite(par_step[[i]])) par_step[[i]] = .Machine$double.xmax
if (par_step[[i]] > (par_hi[[i]] - par_lo[[i]])) par_step[[i]] <- par_hi[[i]] - par_lo[[i]]
}
}
# Reset number of times parameters were accepted
nacp <- rep(0, numpars)
}
# After numpars * ns * nt, reduce temperature
if (cycles %% temp_cycles == 0) {
temp <- temp_red * temp
# NOTE: Goffe et al. restart the search at the current best fit each
# time the temperature drops, but I (CDC) generally don't do this
# Store current maximum lhood in history list
temp_hist_inds <- c(temp_hist_inds, length(lhisthood) + 1)
lhisttemp[length(lhisttemp)+1] <- temp
lhistcycles[length(lhistcycles) + 1] <- cycles / numpars
lhisthood[length(lhisthood) + 1] <- best_lh
temppar <- list()
for (q in 1:length(parnames)) temppar[[parnames[q]]] <- best_par[[parnames[q]]]
lhistpar[[length(lhistpar) + 1]] <- temppar
# Update display
if (show_display) {
likdisplay(lhistcycles, lhisthood, slp, R2, aiccorr, temp, max_iter)
}
# Check to see if likelihood has changed sufficiently to keep going,
# if the user has specified alternate quitting conditions
#if (length(lhisthood) >= min_drops+1) {
if (length(temp_hist_inds) >= min_drops+1) {
if (lhisthood[length(lhisthood)] - lhisthood[temp_hist_inds[length(temp_hist_inds) - min_drops]] < min_change) {
iterate <- FALSE
}
}
}
# Also update display every 100 cycles (but don't update if it just was
# updated)
if (cycles %% display_cycles == 0 && cycles %% temp_cycles != 0) {
lhisttemp[length(lhisttemp)+1] <- temp
lhistcycles[length(lhistcycles) + 1] <- cycles / numpars
lhisthood[length(lhisthood) + 1] <- best_lh
temppar <- list()
for (q in 1:length(parnames)) temppar[[parnames[q]]] <- best_par[[parnames[q]]]
lhistpar[[length(lhistpar) + 1]] <- temppar
if (show_display) {
likdisplay(lhistcycles, lhisthood, slp, R2, aiccorr, temp, max_iter)
}
}
if (using == numpars) using <- 1 else using <- using + 1
if ((cycles / numpars) >= max_iter) iterate <- FALSE
# Here's where we skip to if there was a problem calculating likelihood
# I don't want to update the cycle number yet -
# that way we won't miss any scheduled temperature drops or likelihood updates
}
}
# Calculate support limits
par <- list()
for (i in 1:length(parnames)) {
par[[parnames[i]]] <- best_par[[parnames[i]]]
}
limits<-support_limits(model = model, par = par, var = var, source_data = source_data,
pdf = pdf, par_lo = par_lo_copy, par_hi = par_hi_copy, delta = delta,
slimit = slimit)
upper_limit<-limits$upper_limits
lower_limit<-limits$lower_limits
# Update likelihood history one last time
lhistcycles[length(lhistcycles)] <- cycles / numpars
lhisthood[length(lhisthood)] <- best_lh
temppar <- list()
for (q in 1:length(parnames)) temppar[[parnames[q]]] <- best_par[[parnames[q]]]
lhistpar[[length(lhistpar)]] <- temppar
lhisttemp[length(lhisttemp)] <- temp
if (show_display) {
likdisplay(lhistcycles, lhisthood, slp, R2, aiccorr, temp, max_iter)
}
}, finally={
# Do output
# Update the par_step copy array with the values
# from the "flattened" array
for (i in 1:length(par_step_copy)) {
if (length(par_step_copy[[i]]) == 1)
par_step_copy[[i]] <- par_step[[parnames[[i]]]]
else {
for (j in 1:length(par_step_copy[[i]])) {
par_step_copy[[i]][[j]] <- par_step[[paste(parnames[[i]], j, sep="")]]
}
}
}
# Make par once again just the values that were varied, with the best
# values; we did this when we calculated support intervals, but do it
# again in case the user aborted the run
par <- list()
for (i in 1:length(parnames)) {
par[[parnames[i]]] <- best_par[[parnames[i]]]
}
# Find out what PDF was used
if (identical(pdf, dbeta)) pdfname<-"dbeta"
else if (identical(pdf, dbinom)) pdfname<-"dbinom"
else if (identical(pdf, dcauchy)) pdfname<-"dcauchy"
else if (identical(pdf, dchisq)) pdfname<-"dchisq"
else if (identical(pdf, dnorm)) pdfname<-"dnorm"
else if (identical(pdf, dexp)) pdfname<-"dexp"
else if (identical(pdf, df)) pdfname<-"df"
else if (identical(pdf, dgamma)) pdfname<-"dgamma"
else if (identical(pdf, dgeom)) pdfname<-"dgeom"
else if (identical(pdf, dhyper)) pdfname<-"dhyper"
else if (identical(pdf, dlnorm)) pdfname<-"dlnorm"
else if (identical(pdf, dlogis)) pdfname<-"dlogis"
else if (identical(pdf, dnbinom)) pdfname<-"dnbinom"
else if (identical(pdf, dnorm)) pdfname<-"dnorm"
else if (identical(pdf, dpois)) pdfname<-"dpois"
else if (identical(pdf, dt)) pdfname<-"dt"
else if (identical(pdf, dunif)) pdfname<-"dunif"
else if (identical(pdf, dweibull)) pdfname<-"dweibull"
else if (identical(pdf, dwilcox)) pdfname<-"dwilcox"
else pdfname<-"User-defined function"
# Find out what the user called the model
model_name<-""
base_all<-ls(.GlobalEnv)
for (i in 1:length(base_all)) {
if (identical(get(base_all[i], pos=.GlobalEnv),model)) {
model_name<-base_all[i]
}
}
# Calculate the standard errors, if requested
if (hessian) {
# tryCatch ({
# Flatten par, which has our best parameters
forhess <- numeric()
for (i in 1:length(par)) {
forhess <- c(forhess, par[[i]])
}
HResults <- fdHess(pars=forhess, fun=likeli_4_fdHess, model = model, par = par,
var = var, source_data = source_data, pdf = pdf)
if (!is.na(det(HResults$Hessian)) && det(HResults$Hessian) != 0 &&
kappa(HResults$Hessian) < 1E10) {
var_covar_mat = solve(-1*HResults$Hessian)
std_errs_flat = sqrt(diag(var_covar_mat))
# Fold std_errs_flat back up, if required
std_errs = par
k = 1
for (i in 1:length(par)) {
for (j in 1:length(par[[i]])) {
std_errs[[i]][j] = std_errs_flat[k]
k = k + 1
}
}
} else {
var_covar_mat = "[could not be calculated - Hessian matrix is singular]"
std_errs <- par
for (i in 1:length(par)) {
std_errs[[parnames[i]]] = rep("N/A", times = length(par[[i]]))
}
}
# }, finally={
# var_covar_mat = "[could not be calculated due to errors]"
# std_errs <- par
# for (i in 1:length(par)) {
# std_errs[[parnames[i]]] = rep("N/A", times = length(par[[i]]))
# }
# })
}
else {
var_covar_mat = "[not calculated]"
std_errs <- par
for (i in 1:length(par)) {
std_errs[[parnames[i]]] = rep("N/A", times = length(par[[i]]))
}
}
# To reduce var to manageable size, remove any data frames
if (length(var) > 0)
for (i in 1:length(var)) {if (is.data.frame(var[[i]])) var[[i]]<-"[dataset removed]" }
parhistory <- array(dim=c(length(lhistpar),length(par_lo)))
colnames(parhistory)<-names(par_lo)
for (i in 1:length(lhistpar)) {
parvals<-vector()
for (j in 1:length(lhistpar[[i]])) parvals <- c(parvals,lhistpar[[i]][[j]])
parhistory[i,] <- parvals
}
return(list(best_pars = par, var = var, iterations = cycles / numpars,
source_data = data.frame(source_data, predicted=best_predicted),
par_lo = par_lo_copy, par_hi = par_hi_copy, par_step = par_step_copy,
support_interval_range = slimit,
upper_limits = upper_limit, lower_limits = lower_limit,
initial_temp = initial_temp, temp_red = temp_red, ns = ns, nt = nt,
pdf = pdfname, note = note, model=model_name, std_errs = std_errs,
var_covar_mat = var_covar_mat, max_likeli = best_lh, aic_corr = aiccorr,
aic = (-2.0*best_lh) + (2*numpars), slope = slp, R2 = R2, c = c,
likeli_hist = data.frame(temp = lhisttemp, iter = lhistcycles, likeli = lhisthood, parhistory)))}
)
}
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.