Iwo Jima revisited – how to specify external inputs for dynamic simecol models

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

iwojima1 <- new("odeModel",
    main = function (time, init, parms, ...) {
        x <- init
        p <- parms
        f <- approxTime1(inputs, time, rule = 2)["reinforcements"]
        dame <- f - p["j"] * x["jap"]
        djap <- - p["a"] * x["ame"]
        list(c(dame, djap))},
    parms = c(j=0.0577, a=0.0106),
    times = c(from=0, to=36, by=0.01),
    init = c(ame=0, jap=21500),
    inputs = as.matrix(data.frame(
                day = c(0, 0.999, 1, 1.999, 2, 2.999, 3, 4.999, 
                        5, 5.999, 6, 36),
                reinforcements = c(54000, 54000, 0, 0, 6000, 6000, 0, 0, 
                            13000, 13000, 0, 0))),
    solver = "lsoda") 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

iwojima2 <- new("odeModel",
    main = function (time, init, parms, ...) {
        x <- init
        p <- parms
        dame <- f(time) - p["j"] * x["jap"]
        djap <- - p["a"] * x["ame"]
        list(c(dame, djap)) },
    equations = list(
        f = function(t) {
            if (t<1) return(54000)
            else
              if (t<2) return(0)
              else
                if (t<3) return(6000)
                else
                  if (t<5) return(0)
                  else
                    if (t<6) return(13000)
                    else return(0)}), 
    parms = c(j=0.0577, a=0.0106),
    times = c(from=0, to=36, by=0.01),
    init = c(ame=0, jap=21500),
    solver="lsoda")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

#we read the historical data...
obs <- read.csv2("iwo_jima.csv")
obsv <- which(!(is.na(obs$jap)))
#the weightdf dataframe assigns equal weight to all observations,
#except for the NA values in the japanese column (those get 0 weight)
weightdf <- data.frame(ame=rep(1, nrow(obs)), jap = rep(0, nrow=(obs)))
weightdf[1,"jap"] <- 1
weightdf[nrow(obs),"jap"] <- 1
#and now we fit the model parameters with simecol's fitOdeModel
#including our weightdf as weights makes the PORT routine faster
system.time(
  result_fit1 <- fitOdeModel(iwojima1,
    obstime=obs$time,
    yobs=obs[2:3],
    fn=ssqOdeModel,
    weights=weightdf,
    method="PORT",
    scale=c(1/0.1, 1/0.01),
    lower=c(j=0,a=0))
)
system.time(
  result_fit2 <- fitOdeModel(iwojima2,
    obstime=obs$time,
    yobs=obs[2:3],
    fn=ssqOdeModel,
    weights=weightdf,
    method="PORT",
    scale=c(1/0.1, 1/0.01),
    lower=c(j=0,a=0))
)

Using simecol to simulate a bouncing ball

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

#------------------------------------------------------------------
# Bouncing ball simulation using simecol
# Thomas Petzoldt - 04/09/2010
#------------------------------------------------------------------

library(simecol)
ballode <- new("odeModel",
    main = function (t, y, parms) {
        dy1 <- y[2]
        dy2 <- -9.8
        list(c(dy1, dy2))   },
    parms = 0,
    times = seq(0, 37.5, 0.25),
    init = c(height = 20, v = 0),
    solver = function(t, y, func, parms) {
        root <- function(t, y, parms) y[1]
        event <- function(t, y, parms) {
            y[1] <- 0
            y[2] <- -0.9 * y[2]
            return(y)   }
      lsodar(t, y, func, parms = NULL,events = list(func = event, 
        root = TRUE), rootfun = root)}
)

ballode <- sim(ballode)

# Animation
o <- out(ballode)
png("ball%04d.png", width=300, height=300)
for (i in 1:length(times(ballode)))
    plot(1, o[i, 2], main="", ylab="", xlab="", xaxt="n",
       col="red", cex=2, pch=16, ylim=c(0,20))
graphics.off()
#Convert the *.png image sequence into an animate .gif file
#(Note: requires the ImageMagick command line programs)
system("convert -delay 5 *.png bouncing_ball.gif")

And now, Let's Get Physical5

Notes

1. As a general rule, more elegant and nicer looking code is more efficient in terms of execution time.
2. Overfitting is an important topic in the machine learning context. When we fit a model to observed data, we want the in-sample error to be as small as possible. One way to achieve this is by having more complex models incorporating a greater number of parameters. However, there's a penalty to pay for this, because as we incorporate more parameters into our model, we will almost surely fit the errors that are present in our sample of observations accounting for the peculiarities of that particular sample. In other words, overfitting occurs when we fit the noise in our data sample due to having more parameters in our model and this, while it guarantees a lower in-sample error, will almost surely result in a larger out-of-sample error, which is the error we're really interested in reducing. This is known as the principle of Occam's Razor.
   In our case, we don't have a different number of model parameters in either model and the observed differences in goodness of fit may not even constitute a case of overfitting at all. But if we think about it, the linear interpolation of the data frame mechanism would be more true to life to what occurred in the battle field. In the stepwise function implementation (using the equations slot), you suddenly have the stream of incoming reinforcements cut short to 0 as the simulation transitions from \(t\lt1\) to \(t\geq 1\), whereas in the linear interpolation of the data frame implementation, you have a more gradual reduction in the reinforcement stream. The former, being more true to real-life conditions, is probably why we got a better goodness of fit. I hope you followed me...
3. We will later in this post go into details as to the reasons for modeling the loss of kinetic energy by multiplying the velocity by a factor of 0.9
4. ... Or is almost equal to zero, depending on the error tolerance you define (or predefine in the lsodar function call).
5. Sorry about the cheesy reference to that Olivia Newton-John 1980's song. Anyhow, when looking up the Wikipedia article on her while writing this note, I discovered two interesting facts about her ancestry. On the maternal side, she is descended from Max Born, the famous Nobel Prize atomic physicist. Olivia Newton-John's father was involved in the Enigma project at Bletchley Park, in which the allied mathematicians cracked the Nazi ciphers in World War II. Both facts are tangent on the sort of "physical" that I want to talk about in this section, so let's get physical...
6. The quote that Newton (not to be confused with Newton-John of "Physical"-fame) phrased in English - nanos gigantum humeris insidentes - is actually ascribed to Bernard de Chartres in the 12th Century. As to the invention of Differential Calculus, there's also a view according to which Leibniz invented it independently, but this is too great of a dispute to go into in this footnote.

Bibliography

• Olivia Newton-John. (2016, August 9). In Wikipedia, The Free Encyclopedia. Retrieved August 14, 2016, from https://en.wikipedia.org/wiki/Olivia_Newton-John.
• PETZOLDT, T. and RINKE , K. (2007). "simecol: An Object-Oriented Framework for Ecological Modeling in R". Journal of Statistical Software, 22(9), pp. 1–31. ISSN 1548-7660. http://www.jstatsoft.org/v22/i09.
• R DEVELOPMENT CORE TEAM (2009). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0. http://www.R-project.org.
• Standing on the shoulders of giants. (2016, August 7). In Wikipedia, The Free Encyclopedia. Retrieved August 14, 2016, from https://en.wikipedia.org/wiki/Standing_on_the_shoulders_of_giants.