#############################################	
#	This function plots a histogram for 
#  objects in class samp
#
############################################

hist.samp = function(vals, varname="Variable", adj=2)
{
	x=unclass(vals)
	quartz(varname)
	hist(x,50,freq=F,main=paste("Histogram of", varname), xlab=varname)
	lines(density(x,adjust=adj))
	rug(x)
}
#############################################	
#	This function plots the time history 
#  of a sample in class samp
#
############################################

plot.samp = function(vals,varname="x")
{
	x=unclass(vals)
	quartz(varname)
	plot(x,xlab="Sample Number",ylab=varname)
}

###################################################
#
#	You can "plug in" any of the following 3 proposals
#	Both rq and logdq are needed. The rq and logdq must  
#	match each other on each model (compare the code).
#
###################################################

############## Exact (Gibbs) Proposal #############

rqE = function(model)	{	# Gibbs
	if(model == 1) return(0.5)
	return(rbeta(1,H+1,T+1))
}

logdqE = function(model, param)	{	# Gibbs
	if(model == 1) return(0)
	return(dbeta(param, H+1, T+1, log=TRUE))
}

############## Approximate (Gaussian) Proposal ####

rqA = function(model)	{	# Gaussian
	if(model == 1) return(0.5)
	repeat {
		u = rnorm(1,H/(H+T),sqrt(H*T/(H+T)^3))
		if(u>0 && u<1) return(u)
	}
}

logdqA = function(model, param)	{	# Gaussian
	if(model == 1) return(0)
	mu = H/(H+T)
	sg = sqrt(H*T/(H+T)^3)
	return(dnorm(param, mu, sg, log=TRUE)-		log(pnorm(1.0,mean=mu, sd=sg)-		pnorm(0.0,mean=mu, sd=sg)))
	# Must renormalize the truncated normal distribution
	# or the probabilities will be wrong
	# This is the point of the last two lines.
	# The original code posted was wrong (but the last 
	# two lines evaluate to 1 almost exactly so no error
	# was noticed.
}

##############  Uniform proposal ##################

rqU = function(model)	{	# Uniform
	if(model == 1) return(0.5)
	return(runif(1,0,1))
}

logdqU = function(model, param)	{	# Uniform
	if(model == 1) return(0)
	return(log(dunif(param,0,1)))
}
#
#	Comment: We could have simplified this code to
#	just return the value 0 from this function, but
#	that only works in this particular case, where
#	the proposal is U(0,1). If the proposal had been
#	U(0,2), for example, the log would not evaluate
#	to 0 and we would get different proposals under
#	the two models. The proposal flat(0,1) is not
#	equal to flat(a,b) for arbitrary a and b!!!
#
#	Always check that you can make such simplifications
#	before you make them!
#

############### Log prior ##########################

logdprior = function(model, param)	{
	if(model == 1) return(0)
	return(log(dunif(param,0,1)))

#	Comment: It may seem that the prior on the alternative
#	is flat and can be ignored; but it is actually flat(0,1), 
#	so the log evaluates to 0 only when the prior is U(0,1).  
#	If it were flat(0,2), for example, it would not evaluate  
#	to 0, but to log(0.5), and eliminating the second line  
#	of this calculation would give incorrect results.

}

############### Log beta ###########################
#
#	Note that on a uniform proposal, you can delete the
#	calls to BOTH logdprior AND logdq (if you delete BOTH
#	of them) because IN THIS CASE the prior equals the
#	proposal ON BOTH MODELS and they cancel. Usually this 
#	cannot be done. 
#
#	It is important to be sure that every component
#	of beta_{m|n} is represented, either by actual
#	code or by checking the results of deleting code, 
#	or the sampler will not work properly. I have shown
#	all components; delete code at your peril!
#
#	Calculations are done in logs to protect against
#	underflow
#
####################################################


lbeta = function(model, param, logdq)	{
	x = (
		log(pmodels[model]/qmodels[model])	# p(M)/q(M)
		+ dbinom(H, H+T, param, log=TRUE)	# p(data|param,M)
		+ logdprior(model,param)		# p(param|M)
		- logdq(model, param)			# 1/q(param|M)
	)
		return(x)
}

####################################################
#
#	This function samples M and r in a M-H step
#
####################################################

sample.M.r = function(M,p,proposal)	{
	# Independence sampler: proposal doesn't depend on (M,p)
	Mstar = sample(c(1,2),1,replace=TRUE,qmodels)
	pstar = proposal$rq(Mstar)
	accept = 0
	
	lbeta0 = lbeta(M,p,proposal$logdq)
	lbeta1 = lbeta(Mstar,pstar,proposal$logdq)
	
	u = log(runif(1,0,1))

	if(u < (lbeta1-lbeta0))	{
			M = Mstar
			p = pstar
			accept = 1
	}
	list(M=M,p=p,accept=accept)
}


####################################################
#	
#	This function samples all variables in turn. It
#	first samples M and r in a single M-H step
#	Then we can sample other parameters a, b, ... if
#	they are part of the problem.
#
####################################################

samplen = function(N=1000,proposal)	{
	M = 1
	p = 0.5
	M1 = 0
	models = rep(NaN,N)
	ps = rep(NaN,N)
	
	for(i in 1:N)	{
		sam = sample.M.r(M,p,proposal)
		models[i] = M = sam$M
     ps[i] = p = sam$p
	  if(M==1)	{
		  M1 = M1+1	# Accepts on model 1
	  }
  }
	class(ps)="samp"
	p2=ps[models==2]
	class(p2)="samp"
	cat("Odds on M1 = ", M1/(N-M1),"\n")
	list(model=models,p=ps,p2=p2)
}
	
################ Initializations ####################

#	Data...
H = 60
T = 40

#	Prior on models...
pmodels = c(0.5,0.5)

#	Objects naming the different proposals
UNIFORM = list(rq=rqU,logdq=logdqU)
APPROX = list(rq=rqA, logdq=logdqA)
EXACT = list(rq=rqE, logdq=logdqE)

	
#	Ideally, adjust the qmodels vector to agree with the
#	actual posterior probabilities as sampled. This will
#	give the most efficient sampling. You can do this by
#	running a small sample, observing the acceptances, and
#	then adjusting qmodels appropriately.
qmodels = c(0.5,0.5)

#######################################
#	run program and display results
#######################################
run = function(n=1000,proposal=EXACT){
	z = samplen(n,proposal)
	plot(z$p2)
	hist(z$p2)
	cat("pbar = ",mean(z$p2),"\n")
}

run()