################################################################################## ################################################################################## # An R script to use the Pearson chi-squared method # to fit for the rate of exponential rise of # simulated epidemic data # # Author: Sherry Towers # smtowers@asu.edu # Created: Feb 2nd, 2013 # # Copyright Sherry Towers, 2013 # # 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. ################################################################################## ################################################################################## # first read in the simulated data (which were generated under the assumption # that R0=1.5 and gamma=1/5) ################################################################################## delta_t = 7 R0 = 1.5 gamma = 1/5 rho_true = gamma*(R0-1) vyobs_mat = read.table("simulated_exp_R15_alpha02_nb.txt",header=T) vyobs_mat = as.matrix(vyobs_mat) vyobs_mat_pois = read.table("simulated_exp_R15_pois.txt",header=T) vyobs_mat_pois = as.matrix(vyobs_mat_pois) niter = ncol(vyobs_mat) vt = seq(0,(nrow(vyobs_mat)-1)*7,delta_t) ################################################################################## # for each of the simulated data sets, calculate the pearson-chi squared statistic # for various hypothesized values of rho, then find the value of rho that # minimizes the pearson chi-squared statistic # The estimate of the one standard uncertainty on rho is the range # of values of rho that produce a Pearson chi-squared statistic at most # one greater than the minimum value. ################################################################################## par(mfrow=c(1,1)) ldone = 1 vr = numeric(0) evr = numeric(0) vr_pois = numeric(0) evr_pois = numeric(0) for (iter in 1:niter){ if (iter%%10==0) cat("working on sample",iter,"out of",niter,"\n") for (kter in 1:2){ # 1 is NB smeared simulated data, 2 is Pois smeared simulated data num = 1000 # try num different values of vrho min_rho = 0.0001 vrho = seq(min_rho,(rho_true+0.25),length=num) vpears = numeric(0) for (jter in 1:num){ ydata = vyobs_mat[,iter] if (kter==2) ydata = vyobs_mat_pois[,iter] ypred = (1/vrho[jter])*exp(vrho[jter]*vt)*(1-exp(-vrho[jter]*delta_t)) ypred=ypred[ydata>5] ydata=ydata[ydata>5] ypred = ypred*sum(ydata)/sum(ypred) # It is important to normalize!!! vpears = append(vpears, sum((ydata-ypred)^2/ypred) ) } rho_est = vrho[which.min(vpears)] d=vrho[abs(vpears-min(vpears))<=1.0] d = range(d) delta_rho_est = 0.5*(d[2]-d[1]) if (ldone==0){ xmin = rho_est-2*delta_rho_est xmax = rho_est+2*delta_rho_est vpears = vpears[vrho>xmin&vrhoxmin&vrhomin_rho){ # check to make sure we didn't bump up against the min rho hypothesis if (kter==1){ vr = append(vr,rho_est) evr = append(evr,delta_rho_est) }else{ vr_pois = append(vr_pois,rho_est) evr_pois = append(evr_pois,delta_rho_est) } } } } n=(vr-rho_true)/evr n_pois=(vr_pois-rho_true)/evr_pois if (ldone==2) readline("Press to continue...") par(mfrow=c(2,1)) hist(n,xlab="Normalized residuals from Pearson chi-sq fit to NB smeared simulations",ylab="\043 per bin",main="") hist(n_pois,xlab="Normalized residuals from Pearson chi-sq fit to Pois smeared simulations",ylab="\043 per bin",main="") cat("Mean and standard deviation of normalized residuals from fit to NB smeared simulated data:",mean(n),sd(n),"\n") cat("Mean and standard deviation of normalized residuals from fit to Poisson smeared simulated data:",mean(n_pois),sd(n_pois),"\n")