In jmzobitz/demodelr: Simulating Differential Equations with Data

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Here we are going to explore the likelihood function to try to determine the best set of values for a two parameter dataset using the function plotLikelihood.

First you should open up a blank RStudio document. Be sure to have the library MAT369Code installed.

# Uncomment the next two lines if you need to either update or install the MAT 369 Code
# library(devtools)
# install_github("jmzobitz/MAT369Code",build_vignettes = TRUE)
library(MAT369Code)

Set up

The plotLikelihood function requires three things: Values of your parameters An equation used to compare data (this can be a function or the solution to a differential equation) * An indication of we want to plot the likelihood (the default) or the log likelihood.

We are going to plot the likelihood equation surface for the Gause 1932 dataset, contained in the data frame yeast

First let's make a quick plot of the function:

### Make a quick ggplot of the data
# Contour plot of likelihood function
plotData(yeast,x_label = "Time",y_label = "Volume")

Next we will define the equation used to compare our model in the likelihood. As with euler or systems we have certain lines we can edit:

### This needs to be defined
model <- function(parameters,t){
  with(as.list(c(parameters)), {  # Do not edit this line
    a = log(k/0.45-1) # The value of a is determined by the initial condition
    # In the paper, it sets y(0)=0.45
    yOut = k/(1+exp(a+b*t))
    return(yOut)  
  })
}

Next we will identify the ranges of our parameters we wish to investigate the likelihood: this can be easily varied. The command expand.grid quickly takes all possible combinations of your parameters.

# Identify the ranges of the parameters that we wish to investigate
kParam=seq(5,20,length.out=100)
bParam=seq(-1,0,length.out=100)

# This allows for all the possible combinations of parameters
parameters <- expand.grid(k=kParam,b=bParam)

Then we specify if we want to compute the normal likelihood, or the log-likelihood:

# Define if we want to return the log likelihood function or the regular likelihood

logLikely=TRUE

Now we are ready to visualize the likelihood function! The function plotLikelihood will show the surface and then also report back the value where the likelihood is optimized.

 ### Return the optimum value

plotLikelihood(model,yeast,parameters,logLikely)

Finally you can also make a graph of the optimized parameters against the data:

newParam = c(k=12.72727, b=-0.2424242)
continousTime = seq(0,60,length.out=100)
newVolume <- model(newParam,continousTime)

# 
plotFunction_Data(continousTime,newVolume,yeast,x_label = "Hours",y_label = "Volume")