Hello list:

I am using a Lanchester ODE model of the battle of IwoJima as an excercise in parameter fitting.  In this case, the parameters to estimate are the american and japanese effective rates.  There is day-to-day data of remaining combatant american troops during the 36-day battle.  However, from the Japanese side there is no such data- all that is known is that at day 0 there were approximately 21500 troops and by day 36, for all practical purposes, there were 0 japanese troops still fighting.  What this means is that the corresponding observations in the japanese column are all "NA".  Thus the problem is to fit parameters when some observations are missing.

If i try to fit the parameters by fitOdeModel using the data frame with its NA's, i get an error, as NA's are not numerical values.  By using only the american observations and ommiting the two known japanese observation values, i get a reasonable estimate for the effectiveness rates that seems to be much better than that given in the Engel (1954) paper.  

On the other hand, when I use FME i get parameter estimates that do not lead to a good fit between the model and the historical data.  What am I doing wrong? Could it be that this case is not well suited for FME as the error terms do not meet the assumptions of being normally distributed or homoscedastic?

Bellow i have included the data (which must be put into a file called "iwo_jima.csv" in the working directory), the iwo jima model using fitOdeModel and last, the same model using FME.

Thanks and Regards,
jose romero

-------------------------------------------------------------------------------------------------------------------
DATA (must create file called "iwo_jima.csv" in the working directory
-------------------------------------------------------------------------------------------------------------------
"time";"ame";"jap"
0;0;21500
1;52839;"NA"
2;50945;"NA"
4;56031;"NA"
5;53749;"NA"
6;53155;"NA"
7;65250;"NA"
8;64378;"NA"
9;62874;"NA"
10;62339;"NA"
11;61405;"NA"
12;60667;"NA"
13;59549;"NA"
14;59345;"NA"
15;59081;"NA"
16;58779;"NA"
17;58196;"NA"
18;57259;"NA"
19;56641;"NA"
20;54792;"NA"
21;55308;"NA"
22;54796;"NA"
23;54038;"NA"
24;53938;"NA"
25;53347;"NA"
26;53072;"NA"
27;52804;"NA"
28;52735;"NA"
29;52608;"NA"
30;52507;"NA"
31;52462;"NA"
32;52304;"NA"
33;52155;"NA"
34;52155;"NA"
35;52155;"NA"
36;52140;0

---------------------------------------------------------------------------------------------------------------------
fitOdeModel VERSION
---------------------------------------------------------------------------------------------------------------------
library(simecol)

iwojima <- new("odeModel",
    main = function (time, init, parms, ...) {
        x <- init
        p <- parms
        f <- approxTime1(inputs, time, rule = 2)["refuerzos"]
        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(
                dia = c(0, 0.999, 1, 1.999, 2, 2.999, 3, 4.999, 
                        5, 5.999, 6, 36),
                refuerzos = c(54000, 54000, 0, 0, 6000, 6000, 0, 0, 
                            13000, 13000, 0, 0))),
    solver = "rk4"
) 
#iwojima
iwojima.ajustado <- iwojima

#simular modelo original y registrar resultados
sim1 <- out(sim(iwojima))

#se lee la data observacional para ajustar los parametros del modelo
obs <- read.csv2("iwo_jima.csv")

#ahora ajustamos el modelo
resultado <- fitOdeModel(iwojima,
                         obstime=obs$time,
                         yobs=obs[2],
                         fn=ssqOdeModel,
                         method="PORT",
                         scale=c(1/0.1, 1/0.01),
                         lower=c(j=0,a=0))
resultado

#actualizamos los parametros y volvemos a simular
parms(iwojima.ajustado) <- resultado$par
sim2 <- out(sim(iwojima.ajustado))

#finalmente graficamos
matplot(sim1$time, sim1[2:3], 
        main = "Modelo1 : Engel(1954)",
        xlab = "Día",
        ylab = "Americanos, Japoneses",
        type = "l", lty = c("solid", "solid"), 
        col = c("blue","red"), las = 1)
points(obs$time,obs$ame,col="blue")
matplot(sim2$time, sim2[2:3], 
        main = "Modelo2 : fitOdeModel",
        xlab = "Día",
        ylab = "Americanos, Japoneses",
        type = "l", lty = c("solid", "solid"), 
        col = c("blue","red"), las = 1)
points(obs$time,obs$ame,col="blue")

---------------------------------------------------------------------------------------------------------------------
FME VERSION
---------------------------------------------------------------------------------------------------------------------
library(simecol)

iwojima <- new("odeModel",
    main = function (time, init, parms, ...) {
        x <- init
        p <- parms
        f <- approxTime1(inputs, time, rule = 2)["refuerzos"]
        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(
                dia = c(0, 0.999, 1, 1.999, 2, 2.999, 3, 4.999, 
                        5, 5.999, 6, 36),
                refuerzos = c(54000, 54000, 0, 0, 6000, 6000, 0, 0, 
                            13000, 13000, 0, 0))),
    solver = "rk4"
) 
#iwojima
iwojima.ajustado <- iwojima

#simular modelo original y registrar resultados
sim1 <- out(sim(iwojima))

#se lee la data observacional para ajustar los parametros del modelo
obs <- read.csv2("iwo_jima.csv")

#ahora ajustamos el modelo
library(FME)
fCosto <- function(p) {
    parms(iwojima) <- p
    out <- out(sim(iwojima))
    cost <- modCost(out, obs, weight="std")
    cost
}
times(iwojima) <- obs$time
ajuste <- modFit(f = fCosto, p = parms(iwojima))
summary(ajuste)

#actualizamos los parametros y volvemos a simular
parms(iwojima.ajustado) <- coef(ajuste)
sim2 <- out(sim(iwojima.ajustado))

#finalmente graficamos
matplot(sim1$time, sim1[2:3], 
        main = "Modelo1 : Engel(1954)",
        xlab = "Día",
        ylab = "Americanos, Japoneses",
        type = "l", lty = c("solid", "solid"), 
        col = c("blue","red"), las = 1)
points(obs$time,obs$ame,col="blue")
matplot(sim2$time, sim2[2:3], 
        main = "Modelo2 : FME",
        xlab = "Día",
        ylab = "Americanos, Japoneses",
        type = "l", lty = c("solid", "solid"), 
        col = c("blue","red"), las = 1)
points(obs$time,obs$ame,col="blue")





	[[alternative HTML version deleted]]

