# Code to find the optimal amount of thinning for a Markov sample
# It costs 1 unit to advance the chain
# and theta units to evaluate the function.
# The autocorrelation is rho.

effk = function(k,theta,rho){
# Asymptotic efficiency of thinning factor k vs using k=1
# Compute and exponentiate log( effk )
# NB: log1p( x ) = log( 1+x )
t1 = log1p(theta) - log(k+theta)
t2 = log1p(rho) - log1p(-rho)
t3 = log1p(-rho^k) - log1p(rho^k)

exp( t1 + t2 + t3 )
}

leffkprime = function(k,theta,rho){
# Log of asymptotic efficiency at thinning factor k.
# It ignores terms that do not depend on k.

if( any( rho!=0 ) & any( abs(rho^k) == 0) ){
# Basic detection of underflow while still allowing rho=0
  badk = min( k[abs(rho^k)==0] )
  msg = paste("Underflow for k >=",badk,sep=" ")      
  stop(msg)
}
- log(k+theta) + log1p(-rho^k) - log1p(rho^k)
}

getkmax = function(theta,rho){
# Find an upper bound for the optimal thinning fraction k
if( theta<0 )stop("Negative theta")
if( rho<0 )stop("Negative rho")
if( rho >=1 )stop("rho too close to one")

m=1
while( leffkprime(2*m,theta,rho) > leffkprime(m,theta,rho) )
  m = m*2    

2*m
}

kopt = function(theta,rho,klimit=10^7){
# Find optimal k for the given theta and rho
# Stop if kmax is too large. That usually
# means that theta is very large or rho is very nearly one

kmax = getkmax(theta,rho)    
if( kmax > klimit ){
  msg = paste("Optimal k too expensive. It requires checking",kmax,"values.")
  stop(msg)
}
leffvals = leffkprime( 1:kmax,theta,rho )
which.max(leffvals)

}

kopttable = function( thvals = 10^c(-3:3), rhovals = c(.1,.5,1-10^-c(1:6))){

T = length(thvals)
R = length(rhovals)

bestk = matrix(0,T,R)
row.names(bestk) = thvals
colnames(bestk) = rhovals
effbk = bestk
fracf = bestk

for( i in 1:T )
for( j in 1:R ){
  theta = thvals[i]
  rho   = rhovals[j]
  bestk[i,j] = kopt(theta,rho)
  effbk[i,j] = leffkprime(bestk[i,j],theta,rho)-leffkprime(1,theta,rho)
  effbk[i,j] = exp(effbk[i,j])
}

list( bestk=bestk, effbk=effbk )
}