In jmzobitz/MAT369Code: Simulating Differential Equations with Data

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

In order to make this solver work method work, you will need three things:

• Data: A spreadsheet file of the data you wish to use (or built in data).
• Model: A function that we have for our dynamics (this is $\displaystyle \frac{d\vec{y}}{dt} = f(\vec{y},\vec{\alpha},t)$), or an empirical equation $\vec{y}=f(\vec{x},\vec{\alpha})$.
• Bounds: upper and lower bounds on your parameter values.
• (Initial Conditions): If you have a dynamic model.

We will work you through this step by step, with example code.

Setup and package install

We need to update the package MAT369Code as well as install some new packages:

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

library(deSolve)
library(MAT369Code)

Example 1: A differential equation model

We are going to work through an example that relates the The differential equation we are going to solve is: $$ \displaystyle \frac{dR}{dt}=R\cdot (1-R)-aV $$ $$ \displaystyle \frac{dV}{dt}=b\cdot V \cdot (R-V) $$

This model has two parameters $a$ and $b$, which relate to how the resource is used up as visitors come ($a$) and how as the visitors increase, word of mouth leads to a negative effect of it being too crowded ($b$).

Data

For any data assimilation routine you need to have data that defines the code. This can be from a pre-exisiting dataset already in place, or it can be data that we have defined through importing a file in. For this case we are going to use a pre-defined dataset of the number of resources and visitors to a national park:

knitr::kable(parks)

Additionally because we have a differential equation, we need to specify: - the times that these measurements occur - the initial condition - the independent variable (the cost function needs this to work)

input_data <- parks
time = input_data$time  # the output time vector -must match your data for ODE model

initialCondition = c(resources = 0.995, visitors = 0.00167)

independent_var = "time"

Parameters (including upper and lower bounds)

For all of the parameters in our system we need to define the initial guess and the upper and lower bounds of each parameter. This is down through defining three vectors:

# initial guess of parameters
parameters = c(a = 15, b = 4)

# Lower and upper values of parameters
lower_bound <- c(a=10,b=0)
upper_bound <- c(a=30,b=5)

Run Diagnostics (iterations and burn in percentages)

Two additional things need to be defined: how many iterations we wish to run iterations and the percentage of the first few iterations we will not compute due to the "burn-in" period. This number must be between 0 and 1.

iterations = 1000
burn_percentage = 0.3

If you do not specify these values, they will default to 1500 iterations and a burn-in percentage of 0.2.

Define the model utilized and how you will solve this.

What is missing to make the run is a sense of the model we are using and how we will solve it (via numerical methods to solve ordinary differential equations). This can be quite tricky, but I think we are up to the task here. The method will vary if you have a model that is not an ODE, but it still can be done.

# define model dynamics (used in the ODE function) 
model <- function(time, y, params){
  with(as.list(c(params,y)),{
    dRdt = resources*(1-resources)-a*visitors  ## <-- You may edit here (and add additional lines as necessary) 
    dVdt = b*visitors*(resources-visitors)
    list(c(dRdt,dVdt))
  })
}

# define a function that will solve the model
### This will get the names of the variables


solveModel <- function(pars) {
  return(as.data.frame(  ## Do not edit this line
    ode(initialCondition, time, func = model,parms = pars)  ## <-- You may edit here (and add additional lines as necessary)
  )  ## Do not edit this line
  )
}

That is it! All we need to do is to run our code:

mcmcEstimate(input_data,independent_var,parameters,lower_bound,upper_bound,iterations,burn_percentage) 

Example 2: An empirical model

We are going to work through and example that relates the phosphorous content in algae (denoted by $x$) to the phosphorous content in daphnia (denoted by $y$) The equation we are going to fit is: $$ \displaystyle y = c \cdot x^{1/\theta} $$

This model has two parameters $c$ and $theta$.

Data

For any data assimilation routine you need to have data that defines the code. This can be from a pre-exisiting dataset already in place, or it can be data that we have defined through importing a file in. For this case we are going to use a pre-defined dataset of the number of resources and visitors to a national park:

knitr::kable(phosphorous)

Additionally because we just have a function, we need to specify: - the independent variable (the cost function needs this to work)

input_data <- phosphorous
independent_var = "algae"

Parameters (including upper and lower bounds)

For all of the parameters in our system we need to define the initial guess and the upper and lower bounds of each parameter. This is down through defining three vectors:

parameters = c(c = 1, theta = 5)

# Lower and upper values of parameters
lower_bound <- c(c=0,theta=0)
upper_bound <- c(c=10,theta=30)

Run Diagnostics (iterations and burn in percentages)

Two additional things need to be defined: how many iterations we wish to run iterations and the percentage of the first few iterations we will not compute due to the "burn-in" period. This number must be between 0 and 1.

iterations = 1000
burn_percentage = 0.3

If you do not specify these values, they will default to 1500 iterations and a burn-in percentage of 0.2.

Define the model utilized and how you will solve this.

To define the model it is a slight variation in what we did for a differential equation model. In this case, it is just one equation, and what gets returned (in the command list(algae=algae,daphnia=daphnia) are the columns of a data frame. While it looks a little weird, we do need to name the columns.

# define model equation
model <- function(y, params){
  with(as.list(c(params,y)),{
    daphnia <- c*algae^(1/theta)
    list(algae=algae,daphnia=daphnia)   # For a fitting function we just need this
  })
}


# define a function that will solve the model
### This will get the names of the variables


solveModel <- function(pars) {
  out <-model(input_data,pars)  ## <-- You may edit here (and add additional lines as necessary)
  return(data.frame(out))  ## Do not edit this line

}

That is it! All we need to do is to run our code:

mcmcEstimate(input_data,independent_var,parameters,lower_bound,upper_bound,iterations,burn_percentage) 

jmzobitz/MAT369Code documentation built on March 4, 2024, 12:27 a.m.