# parameters and return values see the wiki # this returns a list containing named pairs: # fit=mls object decribing best fit # conf=the 95% confidence interval

exp_fit_given_logL ← function(times, data, theta0, model, logL) {

# construct a general negative log-likelyhood function
L <- function() {
	# reconstruct the vector of named pairs of 
	# parameters from the list of arguments
	l=length(theta0)
	theta=c()
	vars=names(theta0)
	for (i in 1:l) {
		item=c(x=get(vars[i]))
		names(item) <- vars[i]
		theta <-c(theta, item)
	}
	# compute the mean from the model
	mean <- model(times=times, data=data, theta=theta)
	# compute the negative log-likelihood
	L = -logL(theta, data, mean)
	# if not a number, return a very small likelihood
	if (is.na(L)) return(1e10)
	return(L)
}
# replace the input arguments of L by the list of parameters
formals(L)<-as.list(theta0)

# mle2
# sometimes confint may find a better solution. 
# In this case we need to refit the model with the better
# solution as a starting paoint, because the parameter
# names in the returned better fit can be wrong
fit=mle2(L, method='Nelder-Mead', start=as.list(theta0), control=list(maxit=1e6*length(theta0), ndeps=1e-4, reltol=1e-10))
# fit=mle2(L, method='BFGS', start=as.list(theta0), control=list(maxit=1e5*length(theta0), reltol=1e-10))
# if there is only one parameter to be fitted
# the returned confidence interval is a vector
# change it to a matrix
return(fit)

}

exp_fit ← function(times, data, theta0, model) {

# we start with negative binomial
# if theta0 contains no r, 
if (is.na(theta0['r']))
	theta0 = c(theta0, r=0.0001)
while (TRUE) {
	no_r = is.na(theta0['r']);
	if (!no_r)
	 	p = exp_fit_given_logL(times, data, theta0, model, loglik_negbin)
	# if the fitted shape parameter r is too small (<0.001)
	# then it is approximately poisson 
	# in this case we use poisson
	if (no_r || coef(p)['r'] < 0.001) {
		# we need to remove r from theta
		if (!no_r) {
			old_theta0 = coef(p)
			theta0 = c()
			for (i in 1:length(old_theta0))
				if (names(old_theta0[i]) != 'r')
					theta0 = c(theta0, old_theta0[i])
		}
		p = exp_fit_given_logL(times, data, theta0, model, loglik_poisson)
	}
	conf = confint(p)
	if (is.vector(conf)) conf = matrix(conf, nrow=1)
	if (is.matrix(conf)) break 
	# if we reach here, confint must have found a 
	# better solution. we need to repeat
	theta0 = coef(conf)
	names(theta0) <- names(coef(p))
}
return(list(fit=p, conf=conf, r=coef(p)['r']))

}

#

loglik_negbin ← function(theta, data, mean) {

# for neative bimonial, mean = r p /(1-p)
r = 1 / theta['r'];
# if r < 0, return a tiny likelihood
if (r<0)
	return(-1e10)
sum(dnbinom(x=data, size=r, mu=mean, log=TRUE))

}

#

loglik_poisson ← function(theta, data, mean) {

sum(dpois(x=data, lambda=mean, log=TRUE))

}

#

model_SIR ← function(times, data, theta) {

SIR <- function(t, x, theta) {
	with(as.list(c(x, theta)), {
		DecreaseS = -1 * beta * S * I / N
		DecreaseI = -1 * alpha * I
		dS = DecreaseS
		dI = -1 * DecreaseS + DecreaseI
		return(list(c(dS, dI)))
	})
}

with(as.list(theta), {
	N = S+I
	y = ode(c(S=S, I=I), times, SIR, c(beta=beta, alpha=alpha, N=N), method='ode45')
	return(-diff(c(N,y[, 'S'])))
})

}