################################################################################## ################################################################################## # An R script with utilities to numerically solve a set of ODE's using # the simple time-difference method (Euler's method) # # Author Sherry Towers # smtowers@asu.edu # Created Oct 29 2014 # # Copyright Sherry Towers, 2014 # # This script is not guaranteed to be free of bugs and/or errors. # # This script can be freely used and shared, as long as the author information # and copyright in the header remains intact. ################################################################################## solve_model=function(V,beta,mu,gamma,delt,N){ S = N I = 0 R = 0 vS = numeric(0) vI = numeric(0) vR = numeric(0) vRnew = numeric(0) for (i in 1:length(V)){ Snew = S - delt*(beta*V[i]+mu*I)*S/N Inew = I + delt*(beta*V[i]+mu*I)*S/N - delt*gamma*I Rnew = R + delt*gamma*I vS = append(vS,Snew) vI = append(vI,Inew) vR = append(vR,Rnew) vRnew = append(vRnew,delt*gamma*I) S = Snew I = Inew R = Rnew S[S<0] = 0 I[I<0] = 0 R[R<0] = 0 } return(list(Xpred=vRnew,maxR=max(vR),vS=vS,vI=vI,vR=vR)) } ################################################################################## ################################################################################## # method to fit a cubic spline to a vector, V, and obtain interpolated values # of V separated at times delt (the original vector is assumed to have values # separated by one day) ################################################################################## fit_spline=function(V,delt){ time = seq(1,length(V)) timeb = seq(delt,length(V),by=delt) bv = smooth.spline(time,V,df=length(V)-1,all.knots=T) Vinter = predict(bv,timeb)$y Vinter[Vinter<0] = 0 return(list(Vinter=Vinter,time_inter=timeb)) } ################################################################################## # Define the Negative Binomial probability density function, with parameters # ypred and alpha # The mean of this NB distribution is ypred, and variance is sigma^2=ypred+alpha*ypred^2 ################################################################################## my_dnbinom=function(yobs,ypred,alpha){ if (0){ # use the explicit formulation m = ypred k = 1/alpha a = log(choose(k+yobs-1,yobs)) b = yobs*log(m/(m+k)) c = -k*log(1+m/k) proba = a+b+c proba = exp(proba) } # with transformation of variables, we can use R's inbuilt function for the PDF of the # NB distribution... size = 1/alpha prob = size/(ypred+size) probb = dnbinom(yobs,size=size,prob=prob) return(probb) } ################################################################################## # the negative binomial negative log-likelihood # from: # Handling Overdispersion with Negative Binomial and Generalized Poisson Regression Models # Noriszura Ismail and Abdul Aziz Jemain # www.sherrytowers.com/negative_binomial.pdf # # see: http://books.google.com/books?id=SKUXe_PjtRMC&pg=PA77&lpg=PA77&dq=negative+binomial+reduces+to+poisson&source=bl&ots=lI2iB5qZI0&sig=phdYp_f9EhNtL2BYZTnOFa0ifMw&hl=en&sa=X&ei=hoYSUaCxBuOhyAGRmoHwCA&ved=0CDkQ6AEwAQ#v=onepage&q=negative%20binomial%20reduces%20to%20poisson&f=false # regression analysis of count data, by Colin Cameron ################################################################################## my_lnbinom = function(yobs,ypred,alpha){ # this is the explicit log likelihood, calculated a la Cameron if (0){ mylike = 0.0 for (yi in yobs[yobs>=2]){ r = seq(1,yi-1,1) mylike = mylike+sum(log(1+alpha*r)) } mylike = mylike - sum((yobs+1/alpha)*log(1+alpha*ypred)) mylike = mylike + sum(yobs*log(alpha)) # Ismail and Jemain says negative! mylike = mylike + sum(yobs*log(ypred)) } # or, we can just take the sums of the logs of the probabilities from the PDF... if (1){ aprob = my_dnbinom(yobs,ypred,alpha) mylike = sum(log(aprob[aprob>0])) if (sum(aprob)==0) mylike = 1e20 } return(-mylike) }