# Gibbs sampler (a toy example)
#
# Usage:
# > Gibbs.toy()
# > Gibbs.toy(beta=1, LOOP=1e4, m=20, base=c("A","G","C","T"), pattern=c("C","A","T"))

### preliminaries

count = function(th, pattern){ # count pattern in th
  m = length(th)
  pn = length(pattern)
  if(m < pn) return(0)
  cnt = 0
  for(i in 1:(m-pn+1)){
    if(all(th[i:(i+pn-1)] == pattern)){
      cnt = cnt + 1
    }
  }
  cnt
}

pi.star = function(th, beta, pattern){ # target distribution
  2^(beta * count(th, pattern))
}

h.func = function(th, v, beta, base, pattern){ # "proposal"
  bn = length(base)
  h = numeric(bn)
  for(i in 1:bn){
  	th[v] = base[i]
    h[i] = pi.star(th, beta, pattern)
  }
  h / sum(h)
}

Gibbs.toy = function(beta=1, LOOP=1e4, m=20, base=c("A","G","C","T"), pattern=c("C","A","T"), th0=NULL){
  if(is.null(th0)) th0 = sample(base, m, rep=TRUE)
  th = th0
  stat = numeric(LOOP)
  for(Li in 1:LOOP){ # Gibbs sampling
    v = sample(1:m, 1)
    h = h.func(th, v, beta, base, pattern)
    th[v] = sample(base, 1, prob=h)
    stat[Li] = count(th, pattern)
  }
  cat("estimate =", mean(stat), "\n")
#  ess = LOOP / sum(abs(acf(stat, plot=FALSE)$acf))
  r = acf(stat, plot=FALSE)$acf[2]
  ess = LOOP * (1 - r) / (1 + r)  # suppose AR(1)
  cat("effective sample size =", ess, "\n")
  cat("standard error =", sqrt(var(stat) / ess), "\n")
  par(mfrow=c(1,2))
  plot(stat, type="l", main="sample")
  acf(stat, main="autocorrelation")
}

# set.seed(20170615)
# Gibbs.toy()
# dev.copy2eps(file="Gibbs_toy.eps", height=3, width=6)

### EOF