Nothing
#####################################################
# likeli
# Calculates likelihood.
# Arguments:
# model = model function for which to calculate likelihood.
# Arguments to this function will be provided from par
# and source_data.
# par = list of parameters. If this were simulated annealing,
# these would be the varying parameters. Doesn't make too
# much sense in this context but this makes it match anneal.
# Each par element name matches an argument in a function
# (either model or pdf). All values listed in par must be
# numeric vectors.
# var = list of other variables and data needed by the model
# and pdf functions, any type as needed.
# source_data = data frame with dependent variable and
# associated independent variables
# pdf = probability density function. To calculate
# negative log likelihood, use a function such as
# dnorm that can calculate log probability.
#
# Author: Lora Murphy, Cary Institute of Ecosystem Studies
# murphyl@caryinstitute.org
######################################################
likeli <- function(model, par, var, source_data, pdf, ...) {
##
## Error checking
##
if (!is.function(model)) {
stop("likeli: model is not a function.\n")
}
if (!is.list(par)) {
stop("likeli: par is not a list.\n")
}
if (!is.list(var)) {
stop("likeli: var is not a list.\n")
}
if (!is.data.frame(source_data)) {
stop("likeli: source_data is not a data frame.\n")
}
if (length(source_data) < 1) {
stop("likeli: source_data contains no data.\n")
}
if (length(source_data[[1]]) < 1) {
stop("likeli: source_data contains no data.\n")
}
if (!is.function(pdf)) {
stop("likeli: pdf is not a function.\n")
}
# Make sure all values in par are numeric
for (i in 1:length(par)) {
if (!is.numeric(par[[i]])) {
stop("likeli: All values in par must be numeric.\n")
}
if (!is.vector(par[[i]])) {
stop("likeli: All values in par must be vectors.\n")
}
}
# Here's where we'll put results of all pre-evaluations
eval_results <- NULL
# This is where to put the predicted
predicted <- NULL
# Put together var and par
var <- c(par,var)
# Process the PDF call
datasets<-NULL
datasets[[1]]<-list(value=source_data, varname="source_data")
pdf_call<-analyze_function(pdf, list(value=var, varname="var"), NULL, datasets, ...)
# Process the model call
# Get the list of arguments for model
model_call<-analyze_function(model, list(value=var, varname="var"), NULL, datasets, ...)
##
## Perform the calculations
##
# 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))) {
stop("likeli: the model produced math errors.\nLikelihood cannot be calculated.\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
}
}
}
eval_results <- as.list(eval_results)
# Fourth: sum over the PDF
lh<-sum (do.call(pdf,pdf_call$call))
lh
}
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.