model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    N <- S + I + R
    dS <- s - d*S - beta*S*I/N
    dI <- beta*S*I/N - (d+r)*I
    dR <- r*I - d*R
    return(list(c(dS,dI,dR)))  
  }) 
}  

p <- c(beta=2,r=1,s=1,d=0.01)
s <- c(S=100,I=1,R=0)
run(20)
plane(xmax=100,ymax=20)

p <- c(beta=2,r=1,s=0,d=0)
s <- c(S=10000,I=1,R=0)
data <- read.table("data/niameyA.txt",header=TRUE)
timePlot(data,draw=points,xlab="Time in weeks",ylab="Measles cases",log="y")

f <- lm(log(I)~time, data=subset(data,time<=8))
coef(f)
s["I"]<- exp(coef(f)[1])
p["beta"] <- 1 + coef(f)[2]
free <- c("S","beta","r")
f <- fit(data,free=free,lower=0,ymax=1200,show="I")

f <- fit(data,free=free,lower=0,ymax=1200,show="I",bootstrap=500)
pairs(f$bootstrap)



# make figures
size <- 5 #inch
pdf("sif.pdf",width=size,height=size)
par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i")
p <- c(beta=2,r=1,s=1,d=0.01)
s <- c(S=100,I=1,R=0)
plane(xmax=100,ymax=10)
dev.off()