likelihood.R

################
# Markov Chain #
################
library("markovchain")

change_char = function(x){
  if (x == "a") return (1)
  else if (x == "c") return (2)
  else if (x == "g") return (3)
  return (4)
}


###################################
# Likelihood Ratio/ Cross Entropy #
###################################

loglikhood <- function(data){
  out = log(mean(data == data[1]))
  hold = sapply(data, change_char)
  P = markovchainFit(data)$estimate
  for (index in 2:length(data)){
    out = out + log(P[hold[index - 1], hold[index]])
  }
  out
}
aaaaana/changepoints documentation built on May 26, 2019, 1:35 p.m.