###################################################
### chunk number 1: 
###################################################
meas = read.table('meas.csv', sep=',', header=TRUE)
names(meas)


###################################################
### chunk number 2: 
###################################################
plot(meas$time, meas$London, type='b')


###################################################
### chunk number 3: 
###################################################
cum.reg = smooth.spline(cumsum(meas$B), cumsum(meas$London), df=2.5)

D = predict(cum.reg)$y - cumsum(meas$London) #The residuals

#We visualize the cumulative regression and the reconstructed susceptibles
plot(cumsum(meas$B), cumsum(meas$London), type='l')
lines(cum.reg)
abline(a=0,b=1)


###################################################
### chunk number 4: 
###################################################
ur = predict(cum.reg, deriv=1)$y
summary(ur)


###################################################
### chunk number 5: 
###################################################
Ic = meas$London/ur


###################################################
### chunk number 6: 
###################################################
seas = rep(1:26, 21)[1:545]
lInew = log(Ic[2:546])
lIold = log(Ic[1:545])
Dold = D[1:545]


###################################################
### chunk number 7: 
###################################################
Smean = seq(0.01, 0.2, by=0.001)*3.3E6


###################################################
### chunk number 8: 
###################################################
llik = rep(NA, length(Smean))


###################################################
### chunk number 9: 
###################################################
for(i in 1:length(Smean)){
  lSold = log(Smean[i] + Dold)
  glmfit = glm(lInew ~ -1 +as.factor(seas) + lIold + offset(lSold))
  llik[i] = glmfit$deviance
  }
plot(Smean, llik)

Smean[which(llik==min(llik))]


###################################################
### chunk number 10: 
###################################################
lSold = log(Smean[which(llik==min(llik))] + Dold)
glmfit = glm(lInew ~ -1 +as.factor(seas) + lIold + offset(lSold))
summary(glmfit)


###################################################
### chunk number 11: 
###################################################
Smean[which(llik==min(llik))]


###################################################
### chunk number 12: 
###################################################
glmfit$coef[27]


###################################################
### chunk number 13: 
###################################################
glmfit$coef[1:26]