###########################################
dat<-read.csv("commune.measles.csv",header=F)
par(mfrow=c(1,1))
matplot(dat,type="b",pch=19,xlab="biweek",ylab="cases",
 col=c(1,2,4),lty=c(1,1,1),ylim=c(0,1200))
abline(v=12,lty=2,lwd=2,col=grey(.5))
text(12,1100,labels="MSF campaign",pos=2)
	
###########################################

likelihood<-function(par,vecs){
	S0<-par[1]
	beta<-par[2]
	I<-vecs$I
	S<-floor(S0-cumsum(I[1:(length(I)-1)]))
	p<-1-exp(-beta*(I[1:(length(I)-1)]))	
	L<-dbinom(I[2:length(I)],S,p,log=TRUE)	
	Lik<-sum(L,na.rm=TRUE)
	# note that this returns the log-likelihood
}

###########################################

#first read in the data
dat<-read.csv("commune.measles.csv",header=F)
#choose the time series for the Commune 1
I<-dat[,1]
#make the data list to pass to the likelihood function (above)
vecs<-list(I=I)

#set the number of MCMC iterations
iter<-100000		

#set storage vectors
S0<-numeric(iter)
beta<-numeric(iter)

#choose initial values
S0[1]<-floor(1.9*sum(I))
beta[1]<-1e-4

#set step sizes for proposing new candidate values of parameters
S0.step<-10
beta.step<-5e-6

for(i in 2:iter){
	
	#propose a new value for S0 (note that I've restricted 
	#proposals so that S0 must be greater than the sum of reported cases)
	S0.tmp<-0
	while(S0.tmp<sum(I)){
	  S0.tmp<-S0[i-1]+sample(c(-1,1),1)*rpois(1,S0.step)}
	
	# evaluate the marginal posterior at the candidate value of S0
	# note: these are log likelihoods, so we take the sum, not the product
	numer<- likelihood(c(S0.tmp,beta[i-1]),vecs) + 
	   dunif(S0.tmp,sum(I),10*sum(I),log=T)
	
	#evaluate the marginal posterior at the of the old value of S0
	denom<- likelihood(c(S0[i-1],beta[i-1]),vecs) + 
	   dunif(S0[i-1],sum(I),10*sum(I),log=T)
	
	#Accept or Reject using Metropolis-Hastings
	ifelse(runif(1)<exp(numer-denom),S0[i]<-S0.tmp,
	    S0[i]<-S0[i-1])
	
	#Do the same for beta
	beta.tmp<-0
	while(beta.tmp<1e-9){
	  beta.tmp<-beta[i-1]+sample(c(-1,1),1)*beta.step}
	
	numer<- likelihood(c(S0[i],beta.tmp),vecs) + 
	   dbeta(beta.tmp,1,1,log=T)
	denom<- likelihood(c(S0[i],beta[i-1]),vecs) + 
	   dbeta(beta[i-1],1,1,log=T)
	ifelse(runif(1)<exp(numer-denom),beta[i]<-beta.tmp,
	   beta[i]<-beta[i-1])	
	
	#cat(i,".\n")	#this is a handy counter that will stream 
					#the iteration number on the screen, 
					#so you know that your computer is working
}

###########################################

par(mfrow=c(1,2))
hist(runif(5000,sum(I),10*sum(I)))

hist(rbeta(5000,1,1))

	
###########################################
par(mfrow=c(1,2))
plot(beta,xlab="iteration",main="beta",type="l")
plot(S0,xlab="iteration",main="S0",type="l")
	
###########################################

par(mfrow=c(1,2))
acf(beta[50001:100000],lag.max=1000)
acf(S0[50001:100000],lag.max=1000)
	
###########################################

beta.sub<-beta[seq(50001,100000,by=500)]
S0.sub<-S0[seq(50001,100000,by=500)]

par(mfrow=c(1,2))
hist(beta.sub,main="beta")
abline(v=mean(beta.sub),col=2,lwd=2)
abline(v=c(sort(beta.sub)[c(5,95)]),
   col=4,lwd=2,lty=2)
beta.ci.ts<-c(sort(beta.sub)[c(5,95)])
beta.est.ts<-mean(beta.sub)

hist(S0.sub,main="Initial susceptibles")
abline(v=mean(S0.sub),col=2,lwd=2)	
abline(v=c(sort(S0.sub)[c(5,95)]),
   col=4,lwd=2,lty=2)
S0.ci.ts<-c(sort(S0.sub)[c(5,95)])
S0.est.ts<-mean(S0.sub)

###########################################
# Exercise 1
# for simplicity set the initial values for optim equal to the postior means
# but try additional starting values to prove to yourself that they don't matter
init<-c(floor(mean(S0.sub)),mean(beta.sub))  

# Note: I'm using the argument "control=list(fnscale=-1 . . "
# which tells optim to maximize rather than minimize the likelihood.
mle.fit<-optim(init,likelihood,method="L-BFGS-B",
lower=c(sum(I)+10,1e-4),upper=c(1e4,4e-4),
control=list(fnscale=-1,ndeps=c(1,1e-5)),vecs=vecs)

###########################################
#Exercise 1
init<-c(floor(mean(S0.sub)),mean(beta.sub))
mle.fit<-optim(init,likelihood,method="L-BFGS-B",
lower=c(sum(I)+10,1e-4),upper=c(1e4,4e-4),
control=list(fnscale=-1,ndeps=c(1,1e-5)),vecs=vecs)

layout(matrix(c(1,1,2,3),2,2,byrow=T))

S0.t<-floor(seq(6000,20000,length=40))
beta.t<-seq(.0001,.0004,length=40)

lik<-matrix(NA,40,40)
for(i in 1:40){
for(j in 1:40){
lik[i,j]<-likelihood(c(S0=S0.t[i],
	beta=beta.t[j]),vecs=vecs)
}
}

image(S0.t,beta.t,lik)
contour(S0.t,beta.t,lik,levels=(-200*(2:11)),add=T)
ts.only.est<-c(mean(S0.sub),mean(beta.sub))
points(ts.only.est[1],ts.only.est[2],col=2)
text(ts.only.est[1],ts.only.est[2],
   labels="Bayesian Estimate",pos=1,col=2)
points(mle.fit$par[1],mle.fit$par[2],col=4)
text(mle.fit$par[1],mle.fit$par[2],labels="MLE",
   pos=3,col=4)


hist(beta.sub,main="beta")
abline(v=mle.fit$par[2],col=2,lwd=2)
abline(v=c(sort(beta.sub)[c(5,95)]),col=4,lwd=2,lty=2)


hist(S0.sub,main="S0")
abline(v=mle.fit$par[1],col=2,lwd=2)	
abline(v=c(sort(S0.sub)[c(5,95)]),col=4,lwd=2,lty=2)

###########################################

vac<-read.csv("niamey.vacc.csv",header=T)
# split data into 3 communes
vacc.commune<-split(vac[,3],vac[,2])
# count number of children sampled in each commune
samp.size<-sapply(vacc.commune,length)
# rate of vaccination 
p.vacc<-sapply(vacc.commune,mean,na.rm=T)

###########################################

barplot(p.vacc,
 names.arg=c("Commune 1", "Commune 2", "Commune 3"),
 ylim=c(0,1),ylab="proportion vaccinated")
text(x=1:3,y=p.vacc,labels=as.character(samp.size),pos=3)

###########################################
# Exercise #2
#
#first read in the data
dat<-read.csv("commune.measles.csv",header=F)
#choose the time series for the Commune 1
I<-dat[,1]
#make the data list to pass to the likelihood function (above)
vecs<-list(I=I)

#set the number of MCMC iterations
iter<-100000		

#set storage vectors
S0<-numeric(iter)
beta<-numeric(iter)

#choose initial values
S0[1]<-floor(1.9*sum(I))
beta[1]<-1e-4

#set step sizes for proposing new candidate values of parameters
S0.step<-10
beta.step<-5e-6

for(i in 2:iter){
	
	#propose a new value for S0 (note that I've restricted 
	#proposals so that S0 must be greater than the sum of reported cases)
	S0.tmp<-0
	while(S0.tmp<sum(I)){
	   S0.tmp<-S0[i-1]+sample(c(-1,1),1)*rpois(1,S0.step)}
	
	# evaluate the marginal posterior at the candidate value of S0
	# note: these are log likelihoods, so we take the sum, not the product
	numer<- likelihood(c(S0.tmp,beta[i-1]),vecs) + 
	dpois(S0.tmp,(1-.678)*22432,log=T)
	
	#evaluate the marginal posterior at the of the old value of S0
	denom<- likelihood(c(S0[i-1],beta[i-1]),vecs) + 
	dpois(S0[i-1],(1-.678)*22432,log=T)
	
	#Accept or Reject using Metropolis-Hastings
	ifelse(runif(1)<exp(numer-denom),S0[i]<-S0.tmp,
	S0[i]<-S0[i-1])
	
	#Do the same for beta
	beta.tmp<-0
	while(beta.tmp<1e-9){beta.tmp<-beta[i-1]+
	sample(c(-1,1),1)*beta.step}
	
	numer<- likelihood(c(S0[i],beta.tmp),vecs) + 
	dbeta(beta.tmp,1,1,log=T)
	denom<- likelihood(c(S0[i],beta[i-1]),vecs) + 
	dbeta(beta[i-1],1,1,log=T)
	ifelse(runif(1)<exp(numer-denom),beta[i]<-beta.tmp,
	beta[i]<-beta[i-1])	
	
	#cat(i,".\n")
}
quartz()
par(mfrow=c(2,2))
beta.sub<-beta[seq(50001,100000,by=500)]
S0.sub<-S0[seq(50001,100000,by=500)]

plot(beta.sub,xlab="iteration",main="beta",type="l")
plot(S0.sub,xlab="iteration",main="S0",type="l")


hist(beta.sub,main="beta")
abline(v=mean(beta.sub),col=2,lwd=2)
abline(v=c(sort(beta.sub)[c(5,95)]),
   col=4,lwd=2,lty=2)
beta.ci.vacc1<-c(sort(beta.sub)[c(5,95)])
beta.est.vacc1<-mean(beta.sub)

hist(S0.sub,main="Initial Susceptibles")
abline(v=mean(S0.sub),col=2,lwd=2)	
abline(v=c(sort(S0.sub)[c(5,95)]),col=4,
   lwd=2,lty=2)
S0.ci.vacc1<-c(sort(S0.sub)[c(5,95)])
S0.est.vacc1<-mean(S0.sub)

###########################################

likelihood<-function(par,vecs){
	p.vacc<-par[1]
	beta<-par[2]
	I<-vecs$I
	S0<-vecs$births*(1-p.vacc)
	S<-floor(S0-cumsum(I[1:(length(I)-1)]))
	p<-1-exp(-beta*(I[1:(length(I)-1)]))	
	L<-dbinom(I[2:length(I)],S,p,log=T)	
	Lik<-sum(L,na.rm=TRUE)
}

###########################################

likelihood.vacc<-function(par,vecs){
	p.vacc<-par[1]
	x<-sum(vecs$vacc,na.rm=T)
	N<-sum(!is.na(vecs$vacc))
	dbinom(x,N,p.vacc,log=T)
}

###########################################


#first read in the data
dat<-read.csv("commune.measles.csv",header=F)

#read in the niamey vaccination by commune
vac<-read.csv("niamey.vacc.csv",header=T)
vacc.commune<-split(vac[,3],vac[,2])
samp.size<-sapply(vacc.commune,length)

#select the commune 1 data on cases and vaccination
I<-dat[,1]
vacc<-vacc.commune[[1]]
# note that now "vecs" contains the timeseries 
#AND the vaccine survey data, and the births
vecs<-list(I=I,vacc=vacc,births=22432)  


#set the number of MCMC iterations
iter<-100000		

#set storage vectors
p.vacc<-numeric(iter)
beta<-numeric(iter)

#choose initial values
p.vacc[1]<-.6
beta[1]<-1e-4

#set step sizes for proposing new candidate values of parameters
p.vacc.step<-.001
beta.step<-1e-6

for(i in 2:iter){

	p.vacc.tmp<-p.vacc[i-1]+sample(c(-1,1),1)*p.vacc.step
	numer<- likelihood(c(p.vacc.tmp,beta[i-1]),vecs) + 
	   likelihood.vacc(p.vacc.tmp,vecs) + dbeta(p.vacc.tmp,1,1,log=T)
	denom<- likelihood(c(p.vacc[i-1],beta[i-1]),vecs) + 
	   likelihood.vacc(p.vacc[i-1],vecs) + dbeta(p.vacc[i-1],1,1,log=T)
	ifelse(runif(1)<exp(numer-denom),p.vacc[i]<-p.vacc.tmp,
	p.vacc[i]<-p.vacc[i-1])
	
	beta.tmp<-0
	while(beta.tmp<1e-9){
	    beta.tmp<-beta[i-1]+sample(c(-1,1),1)*beta.step}
	numer<- likelihood(c(p.vacc[i],beta.tmp),vecs) + 
	   dbeta(beta.tmp,1,1,log=T)
	denom<- likelihood(c(p.vacc[i],beta[i-1]),vecs) + 
	   dbeta(beta[i-1],1,1,log=T)
	ifelse(runif(1)<exp(numer-denom),beta[i]<-beta.tmp,
	   beta[i]<-beta[i-1])
	
	#cat(i,".\n")
}

###########################################

par(mfrow=c(1,2))
beta.sub<-beta[seq(50001,100000,by=500)]
p.vacc.sub<-p.vacc[seq(50001,100000,by=500)]
plot(beta.sub,xlab="iteration",main="beta",type="l")
plot(p.vacc.sub,xlab="iteration",
   main="proportion vaccinated",type="l")

###########################################

beta.sub<-beta[seq(50001,100000,by=500)]
p.vacc.sub<-p.vacc[seq(50001,100000,by=500)]

S0.t<-floor(seq(6000,20000,length=40))
beta.t<-seq(.0001,.0004,length=40)

likelihood<-function(par,vecs){
	S0<-par[1]
	beta<-par[2]
	I<-vecs$I
	S<-floor(S0-cumsum(I[1:(length(I)-1)]))
	p<-1-exp(-beta*(I[1:(length(I)-1)]))	
	L<-dbinom(I[2:length(I)],S,p,log=TRUE)	
	Lik<-sum(L,na.rm=TRUE)
	# note that this returns the log-likelihood
}
init<-c(floor(mean(S0.sub)),mean(beta.sub))
mle.fit<-optim(init,likelihood,method="L-BFGS-B",
  lower=c(sum(I)+10,1e-4),upper=c(1e4,4e-4),
  control=list(fnscale=-1,ndeps=c(1,1e-5)),vecs=vecs)

lik<-matrix(NA,40,40)
for(i in 1:40){
for(j in 1:40){
lik[i,j]<-likelihood(c(S0=S0.t[i],beta=beta.t[j]),
   vecs=vecs)
}
}

layout(matrix(c(1,1,2,3),2,2,byrow=T))
image(S0.t,beta.t,lik)
contour(S0.t,beta.t,lik,levels=(-200*(2:11)),add=T)
ts.vacc.est<-c(mean(S0.sub),mean(beta.sub))
points(ts.vacc.est[1],ts.vacc.est[2],col=2)
text(ts.vacc.est[1],ts.vacc.est[2],
    labels="Bayesian Estimate",pos=1,col=2)
points(mle.fit$par[1],mle.fit$par[2],col=4)
text(mle.fit$par[1],mle.fit$par[2],labels="MLE",
   pos=3,col=4)

hist(beta.sub,main="beta")
abline(v=mean(beta.sub),col=2,lwd=2)
abline(v=c(sort(beta.sub)[c(5,95)]),col=4,lwd=2,lty=2)
beta.ci.vacc2<-c(sort(beta.sub)[c(5,95)])
beta.est.vacc2<-mean(beta.sub)

hist((1-p.vacc.sub)*22432,main="Initial Susceptibles")
abline(v=mean(1-p.vacc.sub)*22432,col=2,lwd=2)	
abline(v=c(sort((1-p.vacc.sub)*22432)[c(5,95)]),col=4,
   lwd=2,lty=2)
S0.ci.vacc2<-c(sort((1-p.vacc.sub)*22432)[c(5,95)])
S0.est.vacc2<-mean(1-p.vacc.sub)*22432

###########################################

par(mfrow=c(2,1))
S0.ci<-cbind(S0.ci.ts,
S0.ci.vacc1,
S0.ci.vacc2)

plot(1:3,c(S0.est.ts,S0.est.vacc1,S0.est.vacc2),xlim=c(1,3),ylim=range(S0.ci),
    ylab="S0 confidence intervals",xlab=NA,axes=F,type="b",pch=19,col=2)
axis(2)
axis(1,at=1:3,
   labels=c("ts only","Poisson prior","Vaccine data"))
segments(1:3,S0.ci[1,],1:3,S0.ci[2,])

beta.ci<-cbind(beta.ci.ts,
beta.ci.vacc1,
beta.ci.vacc2)

plot(1:3,c(beta.est.ts,beta.est.vacc1,beta.est.vacc2),xlim=c(1,3),ylim=range(beta.ci),
   ylab="beta confidence intervals",xlab=NA,axes=F,type="b",pch=19,col=2)
   
axis(2)
axis(1,at=1:3,labels=c("ts only","Poisson prior","Vaccine data"))
segments(1:3,beta.ci[1,],1:3,beta.ci[2,])


###########################################
#Exercise 3


#first read in the data
dat<-read.csv("commune.measles.csv",header=F)

#read in the niamey vaccination by commune
vac<-read.csv("niamey.vacc.csv",header=T)
vacc.commune<-split(vac[,3],vac[,2])
samp.size<-sapply(vacc.commune,length)

#select the commune 1 data on cases and vaccination
I<-dat[,2]
vacc<-vacc.commune[[2]]
vecs<-list(I=I,vacc=vacc,births=19936)  


#set the number of MCMC iterations
iter<-100000		

#set storage vectors
p.vacc<-numeric(iter)
beta<-numeric(iter)

#choose initial values
p.vacc[1]<-.6
beta[1]<-1e-4

#set step sizes for proposing new candidate values of parameters
p.vacc.step<-.001
beta.step<-1e-6

for(i in 2:iter){

	p.vacc.tmp<-p.vacc[i-1]+sample(c(-1,1),1)*p.vacc.step
	numer<- likelihood(c(p.vacc.tmp,beta[i-1]),vecs) +  
	    likelihood.vacc(p.vacc.tmp,vecs) + dbeta(p.vacc.tmp,1,1,log=T)
	denom<- likelihood(c(p.vacc[i-1],beta[i-1]),vecs) + 
	    likelihood.vacc(p.vacc[i-1],vecs) + dbeta(p.vacc[i-1],1,1,log=T)
	ifelse(runif(1)<exp(numer-denom),
	    p.vacc[i]<-p.vacc.tmp,p.vacc[i]<-p.vacc[i-1])
	
	beta.tmp<-0
	while(beta.tmp<1e-9){beta.tmp<-beta[i-1]+sample(c(-1,1),1)*beta.step}
	numer<- likelihood(c(p.vacc[i],beta.tmp),vecs) + 
	    dbeta(beta.tmp,1,1,log=T)
	denom<- likelihood(c(p.vacc[i],beta[i-1]),vecs) + 
	    dbeta(beta[i-1],1,1,log=T)
	ifelse(runif(1)<exp(numer-denom),
	    beta[i]<-beta.tmp,beta[i]<-beta[i-1])
	
	#cat(i,".\n")
}

beta.sub<-beta[seq(50001,100000,by=500)]
p.vacc.sub<-p.vacc[seq(50001,100000,by=500)]
hist(beta.sub,main="beta")
abline(v=mean(beta.sub),col=2,lwd=2)
abline(v=c(sort(beta.sub)[c(5,95)]),col=4,lwd=2,lty=2)

hist((1-p.vacc.sub)*19936,main="Initial Susceptibles")
abline(v=mean(1-p.vacc.sub)*19936,col=2,lwd=2)	
abline(v=c(sort((1-p.vacc.sub)*19936)[c(5,95)]),
   col=4,lwd=2,lty=2)


###########################################
likelihood<-function(par,vecs){
	p.vacc<-par[1]
	beta<-par[2]
	I2<-par[3]
	I<-vecs$I
	I[2]<-I2
	S0<-vecs$births*(1-p.vacc)
	S<-floor(S0-cumsum(I[1:(length(I)-1)]))
	p<-1-exp(-beta*(I[1:(length(I)-1)]))	
	L<-dbinom(I[2:length(I)],S,p,log=T)	
	Lik<-sum(L,na.rm=TRUE)
}

###########################################

dat<-read.csv("commune.measles.csv",header=F)
# note that the epidemic in Commune 3 did not start until the 4th biweek
I<-dat[4:16,3]		
vac<-read.csv("niamey.vacc.csv",header=T)
vacc.commune<-split(vac[,3],vac[,2])
vacc<-vacc.commune[[3]]
vecs<-list(I=I,vacc=vacc,births=6708)

iter<-100000
p.vacc<-numeric(iter)
beta<-numeric(iter)
I2<-numeric(iter)

p.vacc[1]<-.5
beta[1]<-1e-4
I2[1]<-1

p.vacc.step<-.001
beta.step<-5e-6
I2.step<-1


for(i in 2:iter){
I2.tmp<-I2[i-1]+sample(c(-1,1),1)*I2.step
numer<- likelihood(c(p.vacc[i-1],beta[i-1],I2.tmp),vecs) + 
    dpois(I2.tmp,5,log=T)
denom<- likelihood(c(p.vacc[i-1],beta[i-1],I2[i-1]),vecs) + 
    dpois(I2[i-1],5,log=T)
ifelse(runif(1)<exp(numer-denom),I2[i]<-I2.tmp,I2[i]<-I2[i-1])

p.vacc.tmp<-p.vacc[i-1]+sample(c(-1,1),1)*p.vacc.step
numer<- likelihood(c(p.vacc.tmp,beta[i-1],I2[i]),vecs) + 
    likelihood.vacc(p.vacc.tmp,vecs) + dbeta(p.vacc.tmp,1,1,log=T)
denom<- likelihood(c(p.vacc[i-1],beta[i-1],I2[i]),vecs) + 
    likelihood.vacc(p.vacc[i-1],vecs) + dbeta(p.vacc[i-1],1,1,log=T)
ifelse(runif(1)<exp(numer-denom),
    p.vacc[i]<-p.vacc.tmp,p.vacc[i]<-p.vacc[i-1])


beta.tmp<-0
while(beta.tmp<1e-9){beta.tmp<-beta[i-1]+sample(c(-1,1),1)*beta.step}
numer<- likelihood(c(p.vacc[i],beta.tmp,I2[i]),vecs) + 
    dbeta(beta.tmp,1,1,log=T)
denom<- likelihood(c(p.vacc[i],beta[i-1],I2[i]),vecs) + 
    dbeta(beta[i-1],1,1,log=T)
ifelse(runif(1)<exp(numer-denom),
    beta[i]<-beta.tmp,beta[i]<-beta[i-1])

#cat(i,".\n")

}

###########################################
#Exercise 4
dat<-read.csv("commune.measles.csv",header=F)
# note that the epidemic in Commune 3 did not start until the 4th biweek
I<-dat[4:16,3]		
vac<-read.csv("niamey.vacc.csv",header=T)
vacc.commune<-split(vac[,3],vac[,2])
vacc<-vacc.commune[[3]]
vecs<-list(I=I,vacc=vacc,births=6708)

iter<-100000
p.vacc<-numeric(iter)
beta<-numeric(iter)
I2<-numeric(iter)

p.vacc[1]<-.5
beta[1]<-1e-4
I2[1]<-1

p.vacc.step<-.001
beta.step<-5e-6
I2.step<-1



for(i in 2:iter){
I2.tmp<-I2[i-1]+sample(c(-1,1),1)*I2.step
numer<- likelihood(c(p.vacc[i-1],beta[i-1],I2.tmp),vecs) + 
    dpois(I2.tmp,5,log=T)
denom<- likelihood(c(p.vacc[i-1],beta[i-1],I2[i-1]),vecs) + 
    dpois(I2[i-1],5,log=T)
ifelse(runif(1)<exp(numer-denom),I2[i]<-I2.tmp,I2[i]<-I2[i-1])

p.vacc.tmp<-p.vacc[i-1]+sample(c(-1,1),1)*p.vacc.step
numer<- likelihood(c(p.vacc.tmp,beta[i-1],I2[i]),vecs) + 
    likelihood.vacc(p.vacc.tmp,vecs) + dbeta(p.vacc.tmp,1,1,log=T)
denom<- likelihood(c(p.vacc[i-1],beta[i-1],I2[i]),vecs) + 
    likelihood.vacc(p.vacc[i-1],vecs) + dbeta(p.vacc[i-1],1,1,log=T)
ifelse(runif(1)<exp(numer-denom),
    p.vacc[i]<-p.vacc.tmp,p.vacc[i]<-p.vacc[i-1])


beta.tmp<-0
while(beta.tmp<1e-9){beta.tmp<-beta[i-1]+sample(c(-1,1),1)*beta.step}
numer<- likelihood(c(p.vacc[i],beta.tmp,I2[i]),vecs) + 
    dbeta(beta.tmp,1,1,log=T)
denom<- likelihood(c(p.vacc[i],beta[i-1],I2[i]),vecs) + 
    dbeta(beta[i-1],1,1,log=T)
ifelse(runif(1)<exp(numer-denom),beta[i]<-beta.tmp,beta[i]<-beta[i-1])

#cat(i,".\n")

}

quartz()
par(mfrow=c(2,2))
beta.sub<-beta[seq(50001,100000,by=500)]
p.vacc.sub<-p.vacc[seq(50001,100000,by=500)]

plot(beta.sub,xlab="iteration",main="beta",type="l")
plot(p.vacc.sub,xlab="iteration",main="proportion vaccinated",type="l")
	

hist(beta.sub,main="beta")
abline(v=mean(beta.sub),col=2,lwd=2)
abline(v=c(sort(beta.sub)[c(5,95)]),col=4,lwd=2,lty=2)

hist((1-p.vacc.sub)*6708,main="Initial Susceptibles")
abline(v=mean(1-p.vacc.sub)*6708,col=2,lwd=2)	
abline(v=c(sort((1-p.vacc.sub)*6708)[c(5,95)]),col=4,lwd=2,lty=2)