#Implementation of ARS algorithm requires several functions to
#be defined.  Code at the end of the file was used for a class
#demo.  Currently, a poorly-chosen mesh results in an algorith
#that works quite poorly, but is interesting pedagogically.  

#Compute slopes for initial grid (mesh). Only called once
gridf=function(x,fx){
nxm1=length(x)-1
yim1=fx(x[1])
a=NULL
b=NULL
for(i in 1:nxm1){
yi=fx(x[i+1])
b[i]=(yi-yim1)/(x[i+1]-x[i])
a[i]=yim1-b[i]*x[i]
yim1=yi
}
out=list(a=a,b=b)
return(out)
}

#Compute sampling probabilities given current mesh, intercepts
# and slopes
omegaf=function(x,a,b){
nx=length(x)
omg=NULL
omg[1]=exp(a[1]+b[1]*x[1])/b[1]
omg[2]=exp(a[2]+b[2]*(x[2]-x[1]))/b[2]
omg[nx]=exp(a[nx-2])*(exp(b[nx-2]*x[nx])-exp(b[nx-2]*x[nx-1]))/b[nx-2]
omg[nx+1]=-exp(a[nx-1]+b[nx-1]*x[nx])/b[nx-1]
for(i in 2:(nx-2)){
icept=(a[i-1]-a[i+1])/(b[i+1]-b[i-1])
p1=exp(a[i-1])*(exp(icept*b[i-1])-exp(b[i-1]*x[i]))/b[i-1]
p2=exp(a[i+1])*(exp(b[i+1]*x[i+1])-exp(icept*b[i+1]))/b[i+1]
omg[i+1]=p1+p2
}
return(omg)
}

#Update x, a, b when a new value is added to the mesh
insertf=function(x,xnew,ponew,fx,a,b){
nx=length(x)
fxn=fx(xnew)
b1=(fxn-fx(x[ponew]))/(xnew-x[ponew])
b2=(fx(x[(ponew+1)])-fxn)/(x[(ponew+1)]-xnew)
a1=fxn-b1*xnew
a2=fxn-b2*xnew
a=c(a[1:ponew],a1,a2,a[(ponew+1):(nx-1)])
b=c(b[1:ponew],b1,b2,b[(ponew+1):(nx-1)])
x=c(x[1:ponew],xnew,x[(ponew+1):nx])
out=list(x=x,a=a,b=b)
return(out)
}

#This is from Example 2.26 in the text; it treats each omega(i) as a uniform rv
#rather than the minimum of a couple truncated exponential rv's. It's much
#easier to implement but probably should be refined.
gnsamp=function(omg,x){
nx=length(x)
dn=sum(omg)
isam=sample(1:(nx-2),1,prob=omg/dn)
xout=x[isam]+runif(1)*(x[isam+1]-x[isam])
out=list(xout=xout,isam=isam,wsam=omg[isam]/dn)
return(out)
}

#a given b
funca=function(alph,bet=.25){
# The scaling issues here make me nervous
x=75:94-mean(75:94)
y=c(3,5,7,9,10,18,6,14,11,9,5,11,15,6,11,17,12,15,8,4)
out=alph*sum(y)-exp(alph)*sum(exp(bet*x))-.1*alph^2/2
return(out)
}

#b given a (Currently not used--don't copy it)
funcb(bet,alph){
x=75:94
y=c(3,5,7,9,10,18,6,14,11,9,5,11,15,6,11,17,12,15,8,4)
out=-exp(alph)*sum(x^2*exp(bet*x))-.2
return(out)
}

#Algorithm A.7
#This grid will turn out to be poorly chosen
igrid=seq(0,2.5,by=.5)
gridab=gridf(igrid,funca)
gridab
igridab=omegaf(igrid,gridab$a,gridab$b)
#Note that igridab puts all its weight on one element
#of the mesh
igridab
igridab/sum(igridab)
igridab=igridab[2:(length(igrid)-1)]
xfinal=NULL
itot=1
# Generate 5 sample points
while(itot <= 5){
sampout=gnsamp(igridab,igrid)
iint=sampout$isam
fnsup=min(exp(gridab$a[(iint-1)]+gridab$b[(iint-1)]*sampout$xout),exp(gridab$a[(iint+1)]+gridab$b[(iint+1)]*sampout$xout))
u=runif(1)
u=u*fnsup
finf=exp(gridab$a[sampout$isam]+gridab$b[sampout$isam]*sampout$xout)
if(u<=finf) {xfinal[itot]=sampout$xout;  itot=itot+1} 
if((u<=funca(sampout$xout)) & (u>finf)){xfinal[itot]=sampout$xout; itot=itot+1; updat=insertf(igrid,sampout$xout,sampout$isam,funca,gridab$a,gridab$b);  
gridab=list(a=updat$a,b=updat$b); 
igrid=updat$x; igridab=omegaf(igrid,gridab$a,gridab$b); 
igridab=igridab[2:(length(igrid)-1)]}
}

#Print the sample of size 5. Note that we only sample from the edges
#of the modal element of the mesh.  This occurs because our 
# mesh is too wide, so the upper bound for almost all choices
# of x in the modal element is too large, and we have little probability
#of accepting points in the interior of the modal element of
# the mesh.  We will try this simulation again with a better
# choice for our initial mesh.
xfinal
igrid
#Try again with a better choice of grid
igrid=seq(1.2,1.5,by=.1)
gridab=gridf(igrid,funca)
gridab
igridab=omegaf(igrid,gridab$a,gridab$b)
#Note that igridab puts all its weight on one element
#of the mesh
igridab
igridab/sum(igridab)
igridab=igridab[2:(length(igrid)-1)]
xfinal=NULL
itot=1